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NUMERICAL  SOLUTIONS  OF  THE  TRIPLE-DECK  EQUATIONS 


FOR  SUPERSONIC  AND  SUBSONIC  FLOW  PAST  A HUMP 


M.  Napolitano,  M.J.  Werle  and  R.T.  Davis 


ABSTRACT 

An  efficient  Alternating  Direction  Implicit  numerical 
technique  is  developed  for  solving  the  Triple-Deck  funda- 
mental problem  for  supersonic  and  subsonic  flow.  Flow 
past  a hump  on  a flat  plate  is  considered  as  a test 
case.  Accuracy  and  reliability  of  the  method  are  positively 
tested  versus  linearized  equation  results  for  a small 
hump  height.  Full  nonlinear  flow  patterns  are  feasible, 
including  mild  separation.  Numerical  difficulties  are 
encountered  around  the  reattachment  point  for  large  hump 
heights . 
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I. 


INTRODUCTION 


In  recent  years  the  problem  of  high  Reynolds  number 
laminar  separation  has  received  a systematic  analysis 
through  the  asymptotic  Triple-Deck  theory  of  Stewartson, 

Sychev,  Neiland  and  Messiter  (Refs.  1-7).  Stewartson  (Ref.  3), 
in  particular,  has  shown  that  in  many  cases  of  practical 
interest,  for  sufficiently  small  separation  causing  dis- 
turbances, the  structure  of  the  separation  region  can  be 
analytically  represented  by  a three  layer  model  which  he 
named  the  "Triple-Deck".  Three  regions  are  identified, 

all  having  a characteristic  longitudinal  length  of  order 
3 

e L;  L is  the  characteristic  length  of  the  geometry 
(usually  the  flat  plate  leading  the  disturbance  region) 
and  £ is  the  small  parameter  in  the  perturbation  expansion, 
related  to  the  large  characteristic  flow  Reynolds  number 

g 

by  the  relation  e = 1/Re.  The  outer  region  (outer  deck) 
is  shown  to  have  a characteristic  length  normal  to  the 

3 

body  surface  also  of  order  e L and  is  characterized  by 
inviscid  irrotational  flow.  The  middle  region  (middle  deck) 

. 4 

has  a normal  characteristic  length  of  order  e L and  is 
basically  a displaced  (Blasius)  boundary  layer  flow, 
inviscid  but  rotational  in  nature,  whose  main  effect  is  to 
passively  communicate  the  pressure  interaction  between  the 
outer  and  inner  decks.  This  last  region  (lower  deck)  is 


1 


1 


a viscous  region  of  characteristic  height  of  order  e“*L, 
for  which  simplified  boundary-layer-like  equations  hold 
and  which  is  produced  by  the  necessity  for  the  flow  field 
to  accommodate  the  viscous  boundary  conditions  at  the  body 
surface.  (Figure  la  gives  a sketch  of  the  Triple-Deck 
structure. ) 

After  nondimensionalization  of  the  governing  equations, 
the  entire  flow  field  around  the  separation  region  is 
determined  to  first  order  in  the  small  perturbation  para- 
meter e , by  solving  the  fundamental  Triple-Deck  problem. 

This  consists  of  the  bottom  deck  boundary-layer-like 
equations  coupled  to  a pressure  interaction  lav;,  relating 
the  viscous  displacement  thickness  growth  and  the  inviscid 
outer  deck  pressure  gradient,  according  to  linear  airfoil 
theory. 

Solutions  to  the  Triple-Deck  equations  are  important  for 
two  reasons:  first  they  provide  rational,  complete  solutions 
for  high  Reynolds  number  flows  where  the  classical  boundary 
layer  approach  fails  (e.g.  in  regions  where  separation  and 
reattachment  occur);  second,  they  are  a reliable  test  base 
for  assessing  more  comprehensive  flow  equation  solvers  (e.g. 
numerical  algorithms  for  the  Navier  Stokes  or  Interacting 
Boundary  Layer  equations;  Ref.  10)  . 

Analytical  solutions  for  the  Triple-Deck  problem  can 
only  be  obtained  for  the  linearized  version  of  these  equations, 
that  apply  when  the  perturbation  producing  geometry  is  small 
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! 
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r 


in  the  bottom  deck  scaling  (see  Refs.  8 and  9 for  example). 

Solution  of  the  general  problem  can  only  be  done  numeri- 
cally due  not  only  to  the  nonlinear  nature  of  the  equations, 
but  also  to  the  fact  that  the  pressure  interaction  law 
makes  them  of  the  boundary  value  type  in  the  longitudinal 
direction.  (This  is  true  even  for  the  case  of  supersonic 
Triple-Deck  flow.) 

To  date,  only  a limited  number  of  yuch  solutions  have 
been  obtained  for  this  general  equation  set.  Solutions  are 
available,  for  example,  for  the  free  interaction  phenomenon 
between  a shock  wave  and  a boundary  layer  flow  (Refs.  2,  12) , 
for  supersonic  flow  past  compression  and  expansion  corners 
(Refs.  13-15)  and  forward  and  backward  facing  steps  (Ref.  15). 

Subsonic  flow  at  the  trailing  edge  of  a finite  flat  plate 
has  also  been  studied  by  several  investigators  independently 
(Refs.  7,  16  and  17) . For  this  case  it  has  been  found 
that  some  of  the  results  of  this  Re  00  asymptotic  theory 
are  surprisingly  accurate  for  moderate  and  even  low  values 
of  the  Reynolds  number  and  are  of  great  practical  importance 
for  calculating  the  viscous  drag  on  a finite  flat  plate 
or  on  a thin  airfoil  with  a sharp  trailing  edge  (see 
Ref.  17  for  example).  Finally,  it  has  been  postulated  that 
the  Triple-Deck  equations  are  valid  around  the  separation 
point  even  in  the  case  of  catastrophic  separation  as  in  the 
flow  past  a circular  cylinder  (Refs.  4,  18).  It  has  just 
recently  been  found  that  numerical  solutions  of  this  set  of 

3 
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equations  are  indeed  possible  in  the  circular  cylinder 
case  and  confirm  the  postulated  structure  of  the  separation 
phenomenon  (Ref.  19) . 

While  definite  advances  are  being  made  in  the  under- 
standing of  interaction  flows  through  use  of  the  Triple-Deck 
concept,  progress  is  being  slowed  by  the  general  ineffi- 
ciency of  most  solution  techniques  available  for  this 
problem  in  the  literature.  Due  to  the  boundary  value  and 
highly  nonlinear  nature  of  the  problem,  they  are  usually 
iterative  procedures  with  two  or  more  nested  iteration 
loops,  which  require  a large  amount  of  computational 
effort  (see  Ref.  16  for  example) . For  the  subsonic  case, 
in  particular,  numerical  convergence  is  extremely  slow  and 
delicate:  severe  underrelaxation  is  required  and  in  some 

cases  an  initialization  reasonably  close  to  the  sought  after 
solution  is  necessary  to  avoid  divergent  behavior.  Here 
the  Alternating  Direction  Implicit  (ADI)  relaxation 
technique  successfully  applied  in  many  "elliptic"  two 
dimensional  flow  problems  (see  Refs.  10  or  11  for  example) 
is  taken  into  consideration.  The  ADI  approach  developed  by 
Werle  and  Vatsa  (Ref.  20) , for  solving  the  supersonic 
Interacting  Boundary  Layer  equations,  is  generalized  here 
to  provide  an  efficient  Triple-Deck  solver  for  both  super- 
sonic and  subsonic  flows.  Results  have  been  obtained  for 
three  flow  configurations:  supersonic  and  subsonic  flow 
past  a parabolic  hump  (Fig.  lb) , and  subsonic  flow  past  a 
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quartic  hump  (Fig.  lc) . Comparisons  have  been  made  in  all 


cases  with  available  solutions  (numerical  nonlinear  or 
analytical  linear)  to  assess  the  reliability  and  effi- 
ciency of  the  solution  method.  Wall  shear  and  pressure 
distribution  profiles  are  presented  for  a wide  range  of 
flow  conditions,  including  separation. 
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II.  THE  FUNDAMENTAL  TRIPLE-DECK  PROBLEM 


I . Governing  Equations 

The  present  work  is  primarily  concerned  with  the 
numerical  solution  of  the  fundamental  Triple-Deck 
equations  for  supersonic  as  well  as  subsonic  flow  con- 
ditions. Therefore,  whereas  such  a set  of  equations 
will  be  discussed  here  in  detail,  only  a very  brief 
outline  of  the  theoretical  work  which  led  to  the  formula- 
tion of  the  Triple-Deck  theory  will  be  presented.  For 
the  reader  interested  in  more  details,  Stewartson  (Ref.  3) 
provides  a review  of  all  the  theoretical  advances  and 
results  of  interest  related  to  "Multi-structured  Boundary 
Layers  on  Flat  Plates  and  Related  Bodies"  up  to  1974. 

This  article  is  the  most  complete  description  to  date  of 
the  asymptotic  Triple-Deck  laminar  separation  theory, 
whose  formulation  is  due  to  Stewartson  himself  more  than 
to  anyone  else.  Also  of  interest  is  Jenson's  historical 
outline  (Ref.  15)  of  the  scientific  struggle  that  led 
from  the  basic  understanding  of  the  interaction  between 
a shock  wave  and  a boundary  layer  to  the  formulation  of  a 
fully  consistent  asymptotic  theory  for  laminar  (supersonic 
flow)  separation.  Reference  (15)  also  reviews  the  practical 
problems  to  which  such  a theory  is  applicable. 

In  the  present  study  a simple  physical  problem  will  be 
briefly  discussed  to  provide  a basic  understanding  of  the 
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asymptotic  Triple-Deck  structure  and  of  the  Triple-Deck 

(fundamental  problem)  equations,  which  will  then  be  given 

without  derivation.  Viscous  flow  over  a flat  plate  is 

considered  (see  Figure  la).  The  flat  plate  has,  at  a 

distance  L from  its  leading  edge,  an  excrescence  or  hump 

of  characteristic  length  and  height  of  order  e^L  and  e5L 

g 

respectively;  e is  the  inverse  of  the  characteristic 

g 

Reynolds  number,  e = v/UrL  = 1/Re,  such  that  e <<  1. 

Since  the  hump  is  "asymptotically"  small  in  physical 
variables,  it  is  reasonable  to  expect  that  at  a large 
enough  distance  ahead  of  and  behind  the  hump  the  flow  field 
will  recover  the  undisturbed  Blasius  boundary  layer  flow 
characteristic  of  a flat  plate  geometry.  Near  the  hump, 
though,  the  disturbance  can  be  large  enough  to  cause  flow 
separation  and  a fine  scale  analysis  is  required  to 
determine  the  important  physical  parameters  such  as  wall 
shear  and  heat  transfer  coefficients.  Figure  lb  shows 
the  small  hump  in  the  bottom  deck  reference  system  in  which 
both  its  length  and  height  hF(x)  are  of  order  one.  It 
has  been  shown  (Refs.  3,  8)  by  rigorous  asymptotic  analysis 
that  to  first  order  in  e the  flow  field  around  the 
disturbance  region  (in  this  particular  case  the  hump)  is 
described  by  the  following  set  of  differential  equations: 


Continuity 


9u  3v 
3x  3y 


0 


(2.1) 


' 


f 


7 


and  Momentum 


u 


9_u 

3x 


+ v 


3u  dp 

3y  dx 


with  boundary  conditions: 

u (x, hF (x) ) = v(x,hF (x) ) = 0 

and 

u(x,y  -*■  «)  -*•  y - 6 ~ h F(x)  , 


(2.2) 


(2 . 3a,b) 


(2.4) 


where  6 is  the  viscous  displacement  thickness.  The  present 
set  of  equations  is  formally  identical  to  the  classical 
boundary  layer  ones,  except  for  boundary  condition  (2.4). 
The  important  difference  is  due  to  the  fact  that  the 
pressure  gradient  is  not  an  imposed  known  function  of  x 
but  "interacts"  with  the  total  displacement  height  D 
according  to  linear  airfoil  theory.  That  is,  for  super- 
sonic flow  conditions  we  have: 


d2P  _ d£ 

,2  dx 

dx 


(2.5) 


whereas  in  case  of  subsonic  flow: 


1 

IT 


(2.6) 


where 

D = 6 + hF (x)  (2.7) 

and  the  integral  in  equation  (2.7)  is  a Cauchy  principal 
value  one.  The  problem  described  by  equations  (2.1)  through 
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(2.6)  is  obviously  not  yet  "closed".  Equations  (3.1)  and 
(3.2)  are  a set  of  parabolic  partial  differential  equations 
and  need  initial  conditions,  given  as 


u ( x -*•  - °°,y)  -*•  y (2.8) 

corresponding  to  the  uniform  shear  flow  constituted  by 
the  "bottom  deck"  undisturbed  Blasius  boundary  layer. 
Eauation  (2.5).  or  (2.6)  equivalently,  accounts  for  the 
boundary  value  nature  of  the  problem  and  requires  the 
following  boundary  conditions: 

6 (x  - - oo)  - 0 (2.9) 

and 

6 (x  -*■  00 ) -*■  0 , (2.10) 

which  again  indicate  the  decaying  nature  of  the  distur- 
bance far  ahead  of  and  behind  the  flow  disturbing  hump. 

2 . Mathematical  Character  of  the  Equations 

There  has  been  controversy  about  the  mathematical 
nature  of  the  Triple-Deck  fundamental  problem,  which  only 
very  recently  has  been  generally  acknowledged  to  be  of 
boundary  value  nature.  In  the  supersonic  case,  in  parti- 
cular, the  interaction  law  equation  (2.5)  can  just  as 

well  be  given  as 
, x 

^ = / p dx  , (2.11) 

dx  i x 

and  a marching  solution  technique,  requiring  no  downstream 
boundary  condition,  would  appear  perfectly  legitimate. 

9 
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And,  in  fact,  the  first  numerical  method,  used  to  solve 
the  free  interaction  fundamental  Triple-Deck  problem, 
was  a marching  procedure  (Ref.  2) . Its  failure  of 
determining  a unique  solution  for  the  reverse  flow  region 
was  already  acknowledged  by  the  authors  themselves  (Ref.  2) 
to  be  caused  by  the  lack  of  the  necessary  accounting  for 
upstream  propagation.  A later  computational  technique, 
applied  to  the  same  free  interaction  problem,  was  instead 
able  to  provide  a unique  reverse  flow  solution  by  pre- 
scribing the  asymptotic  velocity  profile  downstream  and 
marching  backward  toward  the  separation  point  (Ref.  12) . 

The  Triple-Deck  equations  for  subsonic  flow  around  the 
trailing  edge  of  a flat  plate  were  also  numerically  solved 
by  a marching  procedure  (Ref.  16) , and  downstream  recovery 
of  the  Goldstein  wake  solution  was  considered  sufficient 
to  prove  that  the  solution  procedure  was  well  posed. 

Later,  Jenson,  Burggraf  and  Rizzetta  (Ref.  13)  developed 
a numerical  scheme  for  solving  the  Triple-Deck  equations 
for  the  ramp  induced  supersonic  laminar  separation.  They 
properly  accounted  for  the  boundary  value  nature  of  the 
mathematical  problem  and  were  capable  of  determining  a 
unique  solution  for  any  value  of  the  normalized  wedge  angle 
even  in  the  presence  of  very  large  separation  bubbles. 
Further,  the  very  similar  Interacting  Boundary  Layer 
equations  were  solved  for  supersonic  flow  past  a compression 
corner  (Ref.  20) . There  it  was  found  that  branching  could 
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not  be  avoided  and  that  the  weak  interaction  solution  could 
not  be  numerically  recovered  unless  the  boundary  value 
nature  of  the  problem  was  taken  into  account  by  specifying 
appropriate  downstream  boundary  conditions  (see  Ref.  20). 

It  is  therefore  very  clear  that  any  numerical  procedure 
developed  for  efficiently  solving  the  Triple-Deck  equations 
must  properly  account  for  the  boundary  value  nature  of  the 
problem  at  hand  in  order  to  avoid  the  divergent  behavior 
due  to  inappropriate  branching  solutions. 

3.  Prandtl  Transposition  Theorem 

A numerical  solution  of  the  Triple-Deck  problem 
would  encounter  the  difficulty  of  prescribing  the  wall 
boundary  conditions  (2.3)  on  the  body  surface  hF(x), 
which  does  not  necessarily  coincide  with  the  grid  points 
in  the  computational  domain.  It  is  well  known,  though, 

(Ref.  21)  that  the  Prandtl  transposition  theorem  allows 
one  to  transpose  the  longitudinal  coordinate  axis  to  the 
body  surface  for  the  case  of  the  classical  Boundary  Layer 
equations,  which  are  formally  identical  to  equations 
(2.1)  and  (2.2).  Also,  Jenson,  Burggraf  and  Rizzetta 
(Ref.  13) , Rizzetta  (Ref.  14) , and  Jenson  (Ref.  15) , 
have  successfully  applied  this  concept  to  a numerical 
solution  of  the  Triple-Deck  equations  for  supersonic  flow 
past  a compression  ramp.  Here,  following  Jenson's  lead, 
a new  set  of  independent  variables  is  introduced  as: 


J 
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and 


x 


X 


(2.13) 


z = y - hF (x)  , (2.14) 

and  the  normal  velocity  component  is  also  redefined 
according  to 

w=  v-  uhF'(x)  . (2.15) 


Applying  chain  rule  differentiation  to  both  the  u and  v 
variables  as  follows 


3u 

3x 


3u 

3x 


(2.16) 


3u  _ ^u 
3y  3z 

and 


3v  _ 3v 
3y  3z 


F (x) 


(2.17) 


(2.18) 


and  substituting  equations  (2.13)  through  (2.18)  into  the 
continuity  and  momentum  equations  (2.1)  and  (2.2),  it  is 
easily  seen  that  these  equations  remain  unchanged.  Their 
final  form  is  given  below  together  with  their  modified 
boundary  conditions. 


+ I”  = o 
»;  »■ 

(2.19) 

3u  . 3u 

dp  3^u 

~ 2 ' 
dx  3z 

(2.20) 

u(x,0)  = w(x,0) 

= 0 , 

(2.21) 

u(x,z-*-°°)  -*■  z - 

6 

(2.22) 
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The  interaction  laws,  equations  (2.5)  and  (2.6),  remain 
identically  the  same  and  are  also  given  here  for  complete- 
ness as 


and 


dx 


(2.23) 


(2.24) 


4 . Jump  Conditions  for  Supersonic  Flow 

In  the  case  where  the  body  surface  under  consideration, 
hF(x),  has  one  or  more  slope  discontinuities,  according 
to  the  theory  of  characteristics,  it  is  possible  that  the 
flow  variables  also  possess  discontinuous  derivatives  at 
the  corresponding  points.  Moreover,  even  if  the  flow 
variables,  u,v  and  their  derivatives  (for  example,  3u/3x) 
were  continuous  across  the  points  of  discontinuity  for 
f'(x),  the  use  of  the  Prandtl  transposition  theorem  would 
induce  discontinuous  behavior  in  the  normal  velocity 
component  w and  the  derivative  of  the  longitudinal  velocity 
component  3u/3x.  In  fact,  assuming  that  the  point  of 

t 1 , 

discontinuity  in  F (x)  is  at  the  origin  x = 0,  that  is 

f'  (0  + ) t*  f'  (OJ  (2.25) 

and  assuming  that  v and  3u/3x  are  continuous  across  x = 0, 
equations  (2.15)  and  (2.16),  written  immediately  at  the 
left  (0_)  and  at  the  right  (0+)  of  the  origin,  directly 
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provide  the  following  jump  conditions  for  w and  3u/3x: 

w (0+)  + uhF ' (0+)  = w(0_)  + uhF1 (0_)  (2.26) 

and 

^ (0.)  - hF1 (0.)  = ^ (0  ) - hF1 (0  ) . (2.27) 

It  is  therefore  clear  that  a numerical  algorithm,  which 
solves  the  Triple-Deck  equations  for  flow  past  a body  with 
points  of  slope  discontinuity  and  also  makes  use  of  the 
Prandtl  transposition  theorem,  needs  to  properly  account 
for  the  jump  conditions  given  in  equations  (2.26)  and 
(2.27)  . 

In  the  case  of  supersonic  flow  past  a sharp  compression 
corner  it  can  be  proved  that  v and  3u/3x  are  indeed 
continuous  through  the  corner  point  and,  by  simply  assuming 
continuous  pressure  through  the  corner  itself,  the  jump 
conditions  (2.26)  and  (2.27)  can  be  formally  derived  and 
the  pressure  gradient  dp/dx  can  be  proven  to  be  continuous 
(Ref.  15) . The  same  analysis  applies  as  well  to  the  case 
of  supersonic  flow  past  a hump  with  discontinuous  slope 
at  both  the  leading  and  trailing  edges.  In  all  these 
cases,  higher  derivatives  are  discontinuous  even  in  the 
original  x,  y bottom  deck  variables  and  the  corner  points 
are  origins  of  Goldstein  sublayers,  which  eliminate  these 
discontinuities.  See  Reference  (15),  which  also  gives  a 
discussion  of  the  structure  of  the  Goldstein  layer  and 
its  implications  as  to  the  accuracy  of  the  numerical  com- 
putations . 
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For  subsonic  flow  past  a hump  with  discontinuous 
slope  (Ref.  8)  as  well  as  for  the  case  of  flow  at  the 
trailing  edge  of  a flat  plate  (Ref.  16)  (for  which  the 


discontinuity  is  in  the  "wall"  boundary  condition) 

3u/3x  as  well  as  dp/dx  and  d6/dx  are  not  continuous  at 
the  points  of  discontinuity.  It  has  instead  been  shown 
that  in  the  bottom  deck  scaling  3u/3x  °°  at  the  RHS 

(Right  Hand  Side)  of  the  slope  discontinuity  point  in  the 
case  of  the  hump  (Ref.  8)  and  that  dp/dx  and  d6/dx  -*■  00 
at  the  RHS  of  the  trailing  edge  of  the  flat  plate  (Ref.  16) . 
Of  course  these  singularities  are  removed  if  a substructure 
of  the  bottom  deck  is  considered  and  the  Goldstein  sublayer 
is  properly  accounted  for.  The  present  study  makes  no 
account  of  such  a scaling.  Therefore  equations  (2.26) 
and  (2.27)  will  be  properly  accounted  for,  at  any  body 
slope  discontinuity  in  the  case  of  supersonic  flow.  For 
subsonic  flow  conditions,  instead,  a particular  treatment 
will  be  needed  for  the  points  of  discontinuity,  to  avoid 
improper  differencing  through  them.  The  details  will  be 
given  later  in  the  section  on  numerical  methods. 
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III.  NUMERICAL  METHOD  FOR  SUPERSONIC  FLOW 


1 . The  ADI  Relaxation  Technique 

The  numerical  procedure,  developed  here  for  solving 
the  fundamental  Triple-Deck  problem  given  in  equations 
(2.19)  through  (2.24),  is  very  similar  to  the  one  success- 
fully applied  by  Werle  and  Vatsa  (Ref.  20)  to  the  solution 
of  the  supersonic  interacting  boundary  layer  equations  in 
Levy  Lees  variables.  A relaxation  type  time  derivative 

for  the  total  displacement  D is  added  to  the  right  hand 

★ 

side  of  the  momentum  equation  (2.20)  to  give  : 


uu  + wu 
X z 


p + u 
X zz 


(3.1) 


and  a two  sweep  relaxation  technique  proceeds  in  time  from 
the  nth  step  to  (n+l)th  step  in  two  increments  of  At/2. 

Over  the  first  half  time  step  (called  the  * step) , the 
momentum  equation  (3.1)  is  written  as 

(uux  + wuz  - uzz)*  = “ P"  + ft  {D*  ' ^ (3’2a) 


whereas  over  the  second  half  time  step  proceeding  to 
tn+^,  it  becomes 


(uu  + wu 
X z 


* 

u ) = 

zz 


• px+1  + It  <°n+1  - D*> 


(3.2b) 


* Note  that  equation  (3.1)  is  not  the  unsteady  bottom  deck 
momentum  equation  which  would  contain  the  longitudinal 
velocity  time  derivative,  8u/9t,  and  therefore  only 
the  "relaxed"  solution  will  have  physical  meaning.  Also, 
the  ' on  the  x coordinate  is  removed  from  now  onward 
for  convenience. 


r. 
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to  give  a formally  second  order  accurate  solution  in  At. 

Equation  (3.2a)  together  with  the  continuity  equation 

(2.19)  constitutes  a set  of  parabolic  partial  differential 

equations  which  is  solved  by  a standard  marching  technique 

* * * 

to  provide  the  values  of  u , w , D , as  later  described 
in  detail.  (Note  that  at  the  * time  level  the  values 
p^  and  Dn  are  known  functions.)  Proceeding  to  the  n+1 
time  level  by  means  of  equation  (3.2b),  its  left  hand  side 
will  be  known  and  the  solution  will  provide  the  values 
for  the  total  displacement,  Dn+\  and  the  pressure  gradient, 
p^+^.  Since  these  unknowns  are  directly  related  by  the 
interaction  law,  equation  (2.5),  and  are  both  independent 
of  the  normal  coordinate,  z,  the  second  sweep  equation 
(3.2b)  is  first  simplified  using  equation  (3.2a)  to  give 

Px  - A(D  - D ) = px  - A(D  - D ) , (3.3) 

where  A = 2/At  has  been  introduced  for  convenience. 

This,  together  with  the  interaction  law,  equation  (2.5), 
finally  gives: 

„n+l  , ,_n+l  * , ^n  , , * _,n.  ,, 

Dxx  - X<D  - D ) > Dxx  - i(D  - D > . (3.4) 

Note  that  the  unknown  second  derivative,  Dn+1,  is  the 

xx 

term  which  properly  allows  for  the  boundary  value  nature 
of  the  Triple-Deck  equations.  The  finite  difference 
representation  of  equation  (3.4)  will  later  be  seen  to  be 
able  to  directly  accommodate  the  proper  upstream  as  well 
as  downstream  boundary  conditions  on  Dn+'*' . 
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2.  Finite  Difference  * Sweep  Equations  and  Superposition 
Technique 

Figure  I shows  the  computational  grid  used  in  the 
present  numerical  scheme.  Indices  j-1,  j and  j+1  are  used 
to  indicate  the  grid  locations  along  the  normal  direction, 
z,  whereas  1,  2 and  3 are  used  instead  in  the  longitudinal 
x one.  (Subscript  i will  also  be  used  for  the  longitudinal 
direction. ) 


Note  that  the  two  convective  terms  in  equation  (3.6) 

have  been  linearized  or  quasilinearized  around  the  previous 

station  1 and  therefore  equations  (3.5)  and  (3.6)  are 

★ * 

both  linear  in  their  unknowns  u_  . and  w0  . , which  also 

r 3 * / 3 

★ 

linearly  depend  on  the  unknown  total  displacement  D2. 

A simple  superposition  technique  can  be  used  to  solve  the 
set  of  equations  without  iteration.  Two  new  dependent 
variables  are  defined  for  each  of  the  original  ones, 
as  follows: 


(3.7) 

(3.8) 


and  are  then  introduced  in  equation  (3.5)  and  (3.6)  to 
give : 


★ ~ * ★ ~ -k  ~ 


(3.10) 
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For  equations  (3.9)  and  (3.10)  to  be  satisfied,  it 

is  necessary  that  all  the  coefficients  of  the  total 
* 

displacement  D2  vanish,  which  condition  produces  the 
following  set  of  equations  for  u.,  w.: 


w . 
D 


w . 


3 


u . 
3 


(3.11) 


and 


+ 


* 

w . _ 

A • 3 ) u 

4Az  ;uj+l 


= - X + 


4Az 


* * 

(U..  ~ U.  . -)  W. 

1/3+1  1/3-1  3 


(3.12) 


Consequently  the  remaining  terms  of  equations  (3.9)  and 
(3.10)  also  vanish,  that  is: 


W.  = W. 
3 


_ Az  Az  „ . Az  * ,*  , * * 

j-1  Ax  °j  Ax  Uj-1  Ax  (ul,j  Ul,j-1)_  Wl,j  Wl,j-1 


(3.13) 


and 


w 


Az 


1/3' 
2 ' 4Az  1 


9 ui  ■ . w1  . 

U.  1 " + -r^U.  + (-4r  - 

3-i  Az2  Ax  j Az2  4Az 


) 


j+1 


*2 


, _n  , , . n ul , -i  1 , * * 

= 1 D-  + (P  ) - -r-J-  + -r-r-  (U.  .^,-U,  . 


'2  ■ 'rx'  2 T 4A¥  vul,j+l_ul,j-l)wj  * (3-14) 


Equations  (3.11)  through  (3.14)  constitute  two  independent 
sets  of  difference  equations  for  which  appropriate  boundary 
conditions  can  be  easily  obtained  by  splitting  those  for 
equations  (3.9)  and  (3.10)  by  means  of  equations  (3.7)  and 
(3.8)  as  follows: 
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The  wall  no  slip,  zero  injection  conditions 


u2  i = w2  i = 0 (3.15a,b) 

provide  the  following  equations 

UL  = uL  = V = wL  = 0 , (3.16a-d) 

★ 

whereas  the  outer  edge  boundary  condition 

<is  u2>j  * 1 <3-17' 

gives 

(|j  U)j  = 1 (3.18) 

and 

<17“>j'°  <3-19) 


Note  that  the  outer  edge  boundary  condition  (2-4) 
is  first  imposed  at  a finite  distance  z = z^  (j=J) 

and  then  replaced  by  its  derivative,  in  order  to  minimize 
the  error  produced  by  imposing  an  asymptotic  boundary 
condition  at  a finite  distance.  The  analytical 
linearized  solution  shows  that  the  perturbation  shear, 
3u/3z,  decays  at  a slower  rate  than  the  velocity,  u. 
Therefore,  if  the  boundary  condition  is  accurate 
enough  for  the  shear,  it  will  be  definitely  so  for  the 
velocity  profile,  but  not  vice  versa.  Also,  the  present 
superposition  technique  produced  numerical  boundary 
layers  near  the  outer  edge  location  for  u^  and  when 

the  velocity  boundary  condition  was  used.  That  is,  the 
outer  edge  region  was  found  to  be  one  of  extremely  large 
gradients  for  u^  and  Uj  in  order  for  them  to  satisfy 

the  appropriate  boundary  conditions.  Such  behavior 
obviously  produced  large  truncation  errors.  The  use  of 
a shear  outer  boundary  condition  was  instead  found  to 
be  fully  satisfactory. 
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Equations  (3.11)  and  (3.12)  with  their  boundary 

conditions  (3.16b,d)  and  (3.19)  as  well  as  equations  (3.13) 

and  (3.14),  with  boundary  conditions  (3.16a,c)  and  (3.18), 

can  be  numerically  solved  using  the  Davis  coupled  scheme 

(Ref.  22),  thus  providing  the  complete  arrays  IK  , VK  , 

Uj , Wj  at  the  * time  level.  (Appendix  A describes  the 

details  of  the  solution  procedure.)  In  order  to  complete 

the  star  sweep  of  the  present  numerical  procedure  it 

remains  to  determine  the  value  of  the  total  displacement 
* 

D2«  This  can  be  done  by  using  the  velocity  outer  edge 
boundary  condition  (2.22)  which  by  means  of  equation 
(3.7)  can  be  written  as 

Uj  + D*  Uj  = Zj  - (D*  - hF2)  . (3.20) 

* 


Note  that  D2  is  the  only  unknown  in  equation  (3.20)  and 
can  immediately  be  solved  for,  to  give 


* zT  - U_  + hF_ 

D2  = -2 i- 2 


(3.21) 


1 + u. 


* 

Many  alternate  ways  of  evaluating  D2  have  been  investigated. 

In  particular  the  one  obtained  as  described  below  has 

produced  a higher  convergence  rate  in  some  cases,  probably 

because  it  explicitly  accounts  for  the  relaxation  time 

derivat  ve,  3D/3t.  The  momentum  equation  is  written  at 

the  outer  edge  in  its  finite  difference  form  as 
* * 


(z. 


* 

- M 


62~°1 

Ax 


1 * 

2 W1,J 


1 * 

2 W2  , J 


-(Px}2+^52  - 


6£)  = 0 


(3.22) 
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r - 


The  grouping  -r-  (w  +w_  ) is  then  evaluated  by  integrating 

4 J-  / vJ  £,  f U 

the  continuity  equation  from  z = z.  to  z = z by  the 

J.  J 

trapezoidal  rule  and  the  following  final  expression  is 
* * 

obtained  for  D2(<$2+hF2) 

* 

+ z - 6 . J u . +u . . 

^IAt-=ri+jL  =xd" 


+ (px>2  +<IJ-61> 


(61+hF2) 


J 

- I 

i = 2 


* * 

U.+U.  ,-u,  • — u.  . i 

1-1 1>3-1 


(3.23) 


It  must  be  stressed  that  conditions  (3.22)  and  (3.23)  are 
fully  consistent  with  each  other  and,  after  convergence, 
produce,  as  expected,  identically  the  same  results. 

A somewhat  similar  superposition  technique  has  been 
used  by  Sykes  in  his  numerical  solution  of  the  Triple  Deck 
fundamental  problem  for  subsonic  flow  past  a quartic 
hump  (Ref.  23) . 


3.  Solution  of  the  n+1  Time  Sweep  Equation 


The  finite  difference  form  of  the  second  sweep  equation 

(3.4)  uses  central  differences  to  approximate  the  second 

derivatives  and  produces  the  following  tridiagonal 

set  of  linear  algebraic  equations. 

1 nn+1  - f-2-  +X)Dn+1  + -L-  Dn+1 
Di-1  (.  2 +A,Di  + 7 2 Di+1 
Ax  Ax  Ax 


Di-i  - 2DJ  + D"+i 


+ X(D"  - 2Di) 


(3.24) 
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Even  though  the  direct  inversion  of  equation  (3.24)  by  the 
Thomas  algorithm  (Ref.  24)  is  a standard  technique,  it 
will  be  briefly  outlined  here  to  provide  a more  complete 
numerical  treatment  of  the  downstream  boundary  condition. 
Equation  (3.24)  is  of  the  form 


A + B D1?+1  + C = H. 

l-l  l i+l  l 


(3.25) 


and  can  be  directly  inverted  by  the  simple  recursion 
formula 

r.n+1  _ _n+l  , - 

D . = E . D . , + F . 

l l l-l  l 


(3.26) 


where  the  coefficients  Ei  and  F\  are  given  respectively 


as 


and 


E = - A 

1 B + C E. 


i+l 


H.  - C F.^. 
p = _1 l+l 

i B + C E. 


(3.27) 


(3.28) 


'i+l 

It  is  then  clear  that  after  the  coefficients  E^  and  F\ 
are  determined  at  the  last  grid  point,  I,  by  means  of  the 
downstream  boundary  condition  on  Dn+^,  all  the  coefficients 
Ei'  Ei  can  evaluated  through  equations  (3.27)  and  (3.28). 
Finally  the  recursion  relation  (3.26)  gives  the  complete 
array  of  unknowns  once  the  initial  boundary  condition 

for  D*j’+'L  is  provided.  In  the  present  study,  the  asymptotic 
upstream  boundary  condition  6(x  -*■  ~°°)  -*•  0 has  been 
directly  applied  to  Dn+'*'  at  a finite  distance  as 
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I 


L 


D^-1  = 0 


(3.29) 


which  was  numerically  proven  to  be  a sufficiently  good 
approximation.  Several  implementations  of  the  asymptotic 
downstream  boundary  condition,  fi  (x  -*■<»)  -*•  0,  were  instead 
used  in  order  to  minimize  the  negative  effect  of  the  known 
slow  algebraic  decay  for  the  displacement  thickness  6. 

Here  the  case  in  which  F(x)  is  identically  zero  around  the 


location  x^.  where  the  downstream  boundary  condition  was 


imposed  is  considered  in  detail,  the  extension  to  the 
compression  ramp  geometry  where  F(x)  = x (h  = a,  Refs.  13-15) 
being  straightforward. 

The  simplest  way  of  imposing  the  downstream  boundary 
condition  is  of  course  to  impose  at  a finite  , sufficiently 
large  downstream  location  Xj 
.n+1 


D = 0 , 


(3.30) 


which  is  identically  satisfied,  for  any  step  size  Ax, 
by  imposing 


Ei  - fj  - 0 


( 3 . 31a ,b) 


The  simple  observation  that  the  derivative  of  an  alge- 
braically decaying  function  decays  at  a faster  rate  than 
the  function  itself,  obviously  suggests  that  the  error  of 
imposing  an  asymptotic  boundary  condition  at  a finite 
distance  will  be  lower  for  a derivative  condition,  that  is 
Nn+1 


dD1 


dx 


-)  = 0 
I 


(3.32) 
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Equation  (3.32)  can  be  easily  seen  to  be  numerically 
satisfied  to  first  order  accuracy  by  an  analysis  similar 
to  that  of  Appendix  A,  but  simpler,  by  setting 

Ej  = 1 and  Fj  = 0 . (3.33a,b) 

Usually  the  displacement  thickness  <S  has  an  algebraic 
decay  downstream,  that  is 

D(x)  = 6(x)  = Kx  Y (for  x >>  1)  (3.34) 

where  y is  determined  by  asymptotic  analysis  whereas  K 
can  be  left  unknown  for  present  purposes.  Equation  (3.34) 
differentiated  once,  provides  the  following  relation 
between  total  displacement  and  its  derivative 

= - I D (for  x >>  1)  . (3.35) 

dx  x 

The  expression  of  dD/dx  given  by  equation  (3.35)  is  then 
used  in  a Taylor  series  expansion  of  Dn+^(x)  around  the 
final  grid  point  location  x to  give 

D = D + D + 0(Ax2)  (3.36) 

1-1  I Xj  I 

which,  written  in  the  familiar  form 

DI  = 1 + yAx/Xj  DI-1  (3.37) 

finally  provides  the  appropriate  values  for  E^,  given 
respectively  as 
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Et  = 1/(1  + and  FT  = 0 . (3.38a,b) 

1 X j -L 

Note  that  more  accurate  implementation  of  this  boundary 
condition  could  be  straightforwardly  obtained  by  using 
higher  derivatives  of  D(x)  in  the  series  expansion  equa- 
tion (3.36),  if  needed.  Also  note  that  this  way  of 
applying  the  downstream  boundary  condition  at  a finite 
distance  x is  of  higher  order  accuracy  in  an  asymptotic 
sense,  with  respect  to  both  the  previously  discussed  ones. 
This  is  because  it  properly  accounts  for  the  leading  term 
in  the  analytical  asymptotic  expansion  for  6 (x) , equation 
(3.34).  The  use  of  this  boundary  condition  in  the  solution 
of  the  second  sweep  equation  (3.4),  though,  provided 
very  little  advantage  with  respect  to  either  equations 
(3.31)  or  (3.33).  In  fact,  the  numerical  solution  for 
Dn+^  could  not  reach  its  analytical  decay  for  values  of 
X-j.  as  large  as  20,  even  though,  a further  increase  of  Xj 
caused  the  final  solution  to  change  only  immediately 
before  Xj  itself.  At  the  same  time,  the  numerically 
evaluated  coefficient  and  were  seen  to  take  on 
constant  values  over  a wide  range  of  the  independent 
variable  x in  the  downstream  region  and  then  suddenly 
change  to  reach  their  imposed  values  at  the  last  station  x^. 
This  phenomenon  seemed  to  indicate  that  the  finite  dif- 
ference equation  (3.24)  possesses  a somewhat  different 
asymptotic  behavior  than  the  corresponding  differential 
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equation.  Allowing  such  a finite  difference  equation  to 
take  on  its  own  asymptotic  behavior,  as  described  below, 
provided  the  most  consistent  and  effective  downstream 
boundary  condition.  Equations  (3.27)  and  (3.28)  are 
rewritten  here  for  convenience 


E = I_A_ 

i B + C E. 


(3.27) 


- _ Hi  ~ C Fi+1 

Fi  " B + C Ei+1 


(3.28) 


If  E^  and  must  become  constant  for  sufficiently  large 

values  of  i,  they  must  satisfy  equations  (3.27)  and  (3.28) 

with  E . , . = E.  and  F . . . = F..  Equation  (3.27)  thus 
l+l  i 1+1  l 

becomes  a quadratic  algebraic  equation  for  E^ , which,  for 
the  values  of  A,  B and  C given  by  Equation  (3.24),  gives 
the  following  two  roots 


, , XAx2  . XAx2* 2 , 

Ei  = 1 + — 2 — ± \ (1  + ~2 ) - 1 


(3.39) 


The  root  corresponding  to  the  minus  sign  is  chosen  in  order 
to  allow  only  for  a decaying  D^+'L.  After,  this  expression 
(3.39)  is  substituted  into  (3.28)  for  E^  and  the  following 
"asymptotic"  expression  for  F\  is  obtained 


Fi  " X 


“si  +|/d  + “si)2  - i 


(3.40) 


Note  that  can  become  a constant  for  large  x^  values 
only  if  the  right  hand  side  of  equation  (3.40),  namely  , 
also  becomes  a constant.  By  analyzing  the  expression 

Dn  - 2D?  + Dn  . * 

H.  = -iii — + X (D?  - 2D.)  (3.41) 

1 Ax  11 


it  is  clear  that  for  large  values  of  x^  the  second  deri- 
vative of  the  total  displacement  Dn  as  well  as  the 

XX 

• * n 

displacement  D^  and  D^  should  go  to  zero.  Due  to  the 
truncation  error  of  the  numerical  scheme  instead, 
was  seen  to  become  a finite  constant  downstream  and  H. 

l 

also  took  on  a constant  value  for  converged  solutions.  That 

«•!  ^ 1 •/(  « "ft 

is,  for  DV  a = Di  = dV,  Hi  = XD^,  as  clearly  appears  from 
equation  (3.41).  Using  at  the  last  grid  point  the 
values  of  Ej  and  F ^ , given  by  equations  (3.39)  and  (3.40) 
respectively,  was  seen  to  eliminate  the  abrupt  change  in 
the  solution  in  the  neighborhood  of  the  final  station  Xj. 
Enforcing  the  proper  asymptotic  decay  of  the  finite 
difference  equations  also  allowed  the  solution  D^  to  reach 
the  same  level  of  accuracy  with  smaller  values  of  x^.. 
Finally,  the  numerical  inaccuracy  leading  to  constant 
downstream  values  for  the  total  displacement  D^  was  seen 
to  diminish  by  reducing  the  step  size  Ax,  and  these  D^ 
values  linearly  vanished  with  Ax  as  they  should. 
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4 . Windward  Differencing  and  Jump  Conditions 

It  is  well  known  that  numerical  techniques,  using 
backward  difference  representations  for  the  longitudinal 
convective  term  (uu^  in  the  present  case)  in  viscous  flow 
computations,  face  numerical  instability  in  regions  of 
reverse  flow.  This  problem,  which  is  basically  due  to  the 
fact  that  no  upstream  information  is  accounted  for, 
therefore  violating  the  physics  of  the  problem,  is  usually 
removed  by  windward  differencing.  The  following  finite 
difference  form  of  the  convective  term  uu  in  the 


momentum  equation  (3.6)  at  the  * time  level  was  therefore 
used  in  the  present  scheme 


(uu 


ui+|ui'  u2-ui  + yiy 


Ax 


*o  * 
u3  ~u2 
Ax 


(3.42) 


Equation  (3.42)  automatically  provides  a backward  repre- 
sentation for  u^  if  u*  > 0 and  a forward  one  in  case 

* * 
u^  < 0.  In  this  case  the  value  u^  must  be  taken  from 

. . *o 

the  previous  time  level  u^  . Also  note  that  in  equation 

(3.42)  the  j index  has  been  dropped  for  simplicity  and 

★ ~ 

U£  stands  for  either  U or  u in  the  previously  discussed 
superposition  technique.  No  upwinding  of  the  wuz  term 
was  in  general  necessary  as  was  the  case  for  Jenson’s 
numerical  technique  (Ref.  15) . 

For  supersonic  flow  past  a sharp  compression  corner 
or  a parabolic  hump,  with  slope  discontinuities,  the 
numerical  procedure  must  be  able  to  accommodate  the 
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discontinuities  or  jumps  in  the  w and  3u/3x,  due  to  the 
use  of  Prandtl  transposition  theorem. 

In  the  previous  chapter  it  was  mentioned  that  the 
pressure  gradient  p and  therefore  D (see  equation  2.5) 

X XX 

are  continuous  across  the  jump  points.  Therefore  no 

special  treatment  is  required  in  the  second  sweep  equation 

• n+1 

and  a central  finite  difference  representation  of  D 

xx 

is  perfectly  legitimate  at  every  grid  point.  To  properly 
treat  the  * sweep  equations,  the  point  of  slope  dis- 
continuity must  always  coincide  with  a grid  point  in  the 
computational  mesh  (see  figure  II) 


Figure  II.  Slope  Discontinuity 

In  this  way,  while  marching  from  point  1 to  point  2 

this  point  will  be  considered  to  lie  immediately  at  the 

* ★ 

left  of  the  slope  discontinuity,  (2=  2_)  and  u 2 and  w2 
can  be  determined  in  the  usual  manner.  In  marching  from 
point  2 to  3,  instead,  point  2 must  be  considered  to  lie 
immediately  at  the  right  of  the  jump  point  (2  = 2+)  in 
order  for  the  finite  difference  representations,  involving 
points  2 and  3,  to  have  the  expected  order  of  accuracy. 


31 


T 


This  can  be  very  easily  obtained  by  simply  updating  the 

★ 

values  of  w2  by  means  of  the  known  jump  conditions 

equation  (2.26)  as  soon  as  they  have  been  evaluated.  That 

* 

is,  after  the  step  1 to  2 has  been  completed,  the  w2 
are  adjusted  to  their  proper  values  in  order  to  proceed 
to  the  2 to  3 step  as  follows 


w2+  w2- 


u2  h(F2+“  F2-) 


(3.43) 


Equation  (3.43)  eliminates  the  need  of  storing  the  two 

•ft 

different  values  of  w2  and  of  selecting  the  appropriate 
one  for  each  computation  step;  it  is  totally  consistent 
with  Jenson's  jump  condition  treatment  (Ref.  15)  and 
was  numerically  found  to  be  fully  satisfactory.  Note 
that  with  the  present  procedure,  since  u2+  = u2_  and 
only  forward  or  backward  difference  representations 
are  used  for  u^,  any  differencing  through  the  jump  is 
avoided;  therefore,  the  discontinuous  behavior  for  u^ 
is  also  properly  accounted  for. 


IV.  NUMERICAL  METHOD  FOR  SUBSONIC  FLOW 


1 . Second  Sweep  Equation 

The  only  difference  between  the  subsonic  and  supersonic 
Triple-Deck  fundamental  problem  is  the  interaction  law, 
equation  (2.6).  Accordingly,  the  numerical  solution  pro- 
cedure for  the  subsonic  case  will  be  very  similar  to  the 
supersonic  one  and  only  the  major  differences  will  be 
described  here  in  detail.  Equations  (3.2a)  and  (3.2b)  are 
valid  in  the  subsonic  case  too  and,  by  equating  their 
right  hand  side,  the  following  expression  for  p11"^  is 

X 

first  obtained 


n+1 


n 


.n+1 


= p“  + X(D“,J-  - 2D  + Dn)  , 


(4.1) 


which  is  identical  to  equation  (3.3).  Then,  introducing 
the  right  hand  side  of  equation  (4.1)  into  the  subsonic 
interaction  law,  equation  (2.6),  the  final  form  of  the 
second  sweep  equation  is  obtained  as  follows 

* 


n+l  _ 1 1 P?  + x<Dn+1  ' 2D’  + dI*> 

"XX  IT  X-£ 


d 5 


(4.2) 


Note  that  equation  (4.2)  is  the  equivalent  of  equation 
(3.4)  and  this  form  is  the  one  which  properly  allows 
for  the  boundary  value  nature  of  the  problem. 

The  numerical  solution  of  the  * sweep  difference 
equations  proceeds  exactly  as  in  the  supersonic  case. 
(Very  minor  modifications  necessary  to  account  for  the 
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unknown  jump  conditions  will  be  indicated  later.)  The 
numerical  solution  of  the  second  sweep  equation  (4.2)  is 
instead  very  peculiar  and  will  be  discussed  here  in  detail. 
The  second  derivative  is  first  replaced  by  a central  finite 
difference  representation,  and  the  Cauchy  integral  on 
the  right  hand  side  is  broken  into  two  parts  to  give: 


DKi  - 2Df1  + *T-i  i ; 

A 2 TT  £ X-C 

Ax  -00 


1-1  ~n+l 


+ i l » 

TT 


-2D  +D 


n 


dC  • 


(4.3) 


Note  that  p (£)  is  a known  function  at  all  the  grid  points 
inside  the  numerical  integration  region  x^-x^.  The  major 
00  P£dC  XI-1  p£d5 

part  of  (f  is  and  is  numerically 

-OO  X_t’  X1  x“^ 

evaluated  by  the  methods  described  in  Appendix  B,  whereas 
the  residuals 
x, 


I1  Vi 

1 X-i 


and  / 
x 


1-1 


x=T~ 


can  be  assessed  if  the  asymptotic  decay  of  p(x)  is  known 

far  ahead  of  and  far  behind  the  region  of  disturbance 

(see  Appendix  D for  details) . Also  note  that  in  equation 

(4.3)  the  Cauchy  integral  of  the  relaxation  time  dependent 

term  has  been  approximated  with  its  value  over  the  region 

* 

of  numerical  integration  only  . This  , because  its  purpose 


* According  to  the  Cauchy  integral  numerical  evaluation  given 
in  Appendix  A,  a Cauchy  integral  can_be  evaluated  either 
at  a grid  point  x^  or  at  a midpoint  x^  = x-^+Ax^  The 
corresponding  regions  of  integration  are  x^  - x^ and 
x^  - Xj  respectively. 
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is  to  allow  the  solution  to  relax  to  its  steady  state  value 
and,  after  this  has  been  reached,  the  integral  vanishes. 
Unfortunately  though,  equation  (4.3)  cannot  be  reduced  to 
a tridiagonal  set  of  linear  algebraic  equations  and  in- 
verted by  Thomas  Algorithm  because  of  the  presence  of  the 
unknown  Dn+^  inside  the  Cauchy  integral.  An  iterative 
procedure  is  then  necessary  to  solve  the  second  sweep  of 
the  two  step  ADI  method,  equation  (4.3).  The  most 
efficient  solution  procedure  which  we  have  been  able  to 
develop  consists  of  extracting  from  the  Cauchy  integral 
the  upper  and  lower  diagonal  entries  for  the  D?+^  array, 

replacing  all  the  others  with  their  previous  iteration 
n+1 

level  values  Dp_.  , and  inverting  the  resulting  tndiagonal 
matrix  by  the  Thomas  Algorithm.  In  more  detail,  the 


integral 


Cl 


<V  - 7 


1-1  „n+l 


-2D  + D 


n 


Xi-S 


(4.4) 


appearing  in  the  right  hand  side  of  equation  (4.3)  is 

written  in  finite  difference  form  according  to  equation 

(B2 ) of  Appendix  B at  any  location  and  the  tridiagonal 

elements  contributions,  namely  ^ £n3  D*?*^,  0,  ~ £n3  D?+^ , 

are  brought  to  the  left  hand  side  of  equation  (4.3).  By 

also  replacing  all  the  other  D^+^  contributors  to  the 

integral  CI(x.)  of  equation  (4.4)  with  their  known  previous 
1 n+1 

iteration  values  Dp  . , equation  (4.3)  is  reduced  to  the 
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following  tridiagonal  finite  difference  equation 


X . _n+l 
- — £n  3)  D . , 
tt  i-l 


Dn+1  + (—•  + 2,n3)Dn^  = H. 

Ax2  1 Ax2  ^ 1+1 


(4.5) 


in  which  the  known  right  hand  side  IT  also  accounts  for 


XI-1  p?  dC 


j.  j.  r v. 

the  A -7T1 


contributor. 


The  inversion  procedure  for  equation  (4.5)  is  now 
identical  to  the  one  used  for  the  second  sweep  equation 

(3.4)  in  the  supersonic  case,  but  it  must  be  repeated 

until  some  convergence  criterion  for  Dn+"''  and  Dn+'*'  is  met. 

1 pi 

As  in  many  previous  studies  (Refs.  16  and  17) , under- 
relaxation was  found  necessary  to  avoid  divergent  behavior. 
At  the  end  of  each  iteration  the  new  value  of 
obtained  by  inverting  the  tridiagonal  matrix,  equation 

(4.5) ,  was  underrelaxed  according  to  the  following 
relation 


Dn+1  *■  0.2  Dn+1  + (1  - 0.2)  Dn+1 
l l pi 

in  which  the  numerical  value  of  the  underrelaxation 


(4.6) 


coefficient  0.2  was  taken  as  in  Ref.  (16).  Other  under- 
relaxation coefficients  were  found  to  provide  convergent 
behavior  for  the  second  sweep  equation,  but  no  systematic 
numerical  study  was  made  in  order  to  find  the  optimal 
value.  It  was  found  that  the  relaxation  parameter  had  to 


be  lower  than  0.5  in  order  to  obtain  convergent  behavior. 


r 


I 

f 


Since  the  present  relaxation  technique  does  not 
provide  a physically  meaningful  transient,  a complete 
convergence  in  the  second  sweep  iterative  procedure  was 
not  necessary  at  every  time  step;  a limit  of  10  was 
therefore  imposed  for  the  number  of  iterations  performed 
at  each  n+1  time  level. 

All  the  treatments  of  the  boundary  conditions  described 
for  the  supersonic  second  sweep  equation  (3.4)  apply  as 
well  to  equation  (4.5).  Allowing  the  coefficients  , 

F\  of  the  recursion  formula,  equation  (3.26),  to  take  on 
their  "asymptotic"  values  was  again  found  to  be  the  most 
satisfactory  approach.  and  F could  then  be  determined 

exactly,  as  in  the  supersonic  case  of  equations  (3.39) 
and  (3.40),  to  be 


E = (1  - Mill  Ax2 ) / ( 1 + MUi  Ax2) 

I IT  IT 


and 


Fi  ’ 


H Ax 


Ejd  + Ax2  - 2) 


(4.7) 


(4.8) 


Fj  was  found  to  become  practically  zero  in  all  converged 
numerical  computations,  but  sometimes  it  grew  without  bound, 
thus  inducing  divergent  behavior  in  the  numerical  procedure. 
Setting  Fj  identically  equal  to  zero  (which  it  is  to  second 
order  accuracy  in  Ax)  was  sufficient  to  remove  this  problem 
without  deteriorating  the  accuracy  of  the  solution. 
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2 . * Sweep  Slope  Discontinuities  Treatment 

For  subsonic  flow  past  a parabolic  hump  with  slope 
discontinuities,  the  linearized  Triple-Deck  equations  (Ref.  8) 
produce  pressure  and  wall  shear  distributions  with  dis- 
continuous gradients  at  any  point  of  slope  discontinuity 
(x  = 0) , that  is 


and 


(§E)  ft  (£E) 

',dx;  r vdx;  ' 

QX  0+  0- 


8x  9t 

(_J£)  u (_W) 

3x  o+  9x 


(4.9) 


(4.10) 


where 


, 3u. 

Tw  3z  z=o 


(4.11) 


Furthermore,  the  wall  shear  gradient  becomes  singular  at 
the  right  hand  side  of  the  slope  discontinuity  point  (Ref.  8) , 
that  is 

3t 

- « . (4.12) 


(— ) 
l3x  ' 


0+ 


It  is  obvious  that,  in  the  absence  of  analytically  prescribed 
jump  conditions,  a special  treatment  of  the  finite  difference 
equations  is  required  to  properly  account  for  such  behavior 
of  the  flow  variables.  The  * sweep  equations  were  treated 
as  follows.  The  jump  point  was  again  made  to  always  coincide 
with  a grid  point  x ^ . In  marching  from  point  x^  to  x^ , x^ 
was  considered  to  lie  immediately  ahead  of  the  jump  and  the 
usual  finite  difference  equations  (3.5)  and  (3.6)  were 
perfectly  legitimate.  In  marching  from  point  X£  to  x^. 
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*2  was  considered  to  lie  immediately  at  the  right  of  the 
jump  point  and  the  following  finite  difference  representa- 
tions were  used  for  the  continuity  equation  and  for  the 


convective  term  wuz  in  the  momentum  equation  respectively. 


* * 
• ~u. 


* * 

w.  .-w. 


1U.  • Ua  1 U A • - ""Ua  1 VY.  . W—  • a 

. Aj  j 1 + i 3 > 2,1-1  _3  yj 3 , J-l  g /.  1 3) 

2 Ax  2 Ax  Az  U 


and 


* * 

* u_  . -u_  . . 

w 2,j-l  2,3-1 

3,j  2Az 


(4.14) 


(Note  that  in  equations  (4.13)  and  (4.14)  the  unknown 

variables  to  which  the  superposition  technique  applies 
★ * 

are  u,  . and  w-  ..)  The  finite  difference  representations 
j , 3 j , j 

(4.13)  and  (4.14)  eliminate  the  presence  of  the  unknown 
* 

w_  . array  and  avoids  differencing  across  the  jump  point. 

* t 3 

They  are  therefore  perfectly  legitimate  as  long  as 

* * 

U2+  - u2_'  that  is  the  u-velocity  is  continuous  across 
the  jump,  which  it  is.  The  modifications  required  for  the 
coefficients  of  the  inversion  algorithm  of  Appendix  A, 
in  order  to  use  equations  (4.13)  and  (4.14)  at  each  grid 
following  a jump  point  are  very  straightforward. 


3 . n+1  Sweep  Slope  Discontinuities  Treatment 

It  has  been  mentioned  that  the  pressure  gradient  px 
is  discontinuous  across  the  points  of  slope  discontinuity. 
The  second  derivative  of  the  total  displacement,  given  as 


D 


xx 


X-S 


(2.6) 
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becomes  singular  in  correspondence  to  the  points  of  dis- 
continuity for  (See  Appendix  E for  the  analytical 

proof.)  It  appears  therefore  improper  to  express  D 
with  a central  difference  at  the  corner  points,  since  the 
Taylor  series  approach  on  which  this  representation  is 
based  evidently  fails.  In  order  to  remove  this  problem 
the  following  approximation  has  been  used  to  model  dp/dx, 
which  provides  a correct  evaluation  of  D in  the  limit 

XX 

of  zero  step  size  and  an  "approximate"  regular  one  for 
finite  values  of  Ax. 

dp/dx 


Figure  III.  Model  for  the  Pressure  Gradient 
Discontinuity . 

Figure  III  shows  an  analytical  function  dp/dx  discontinuous 
across  x^.  In  the  present  numerical  procedure  dp/dx  is 
only  known  at  each  grid  point  x^.  However,  in  order  to 
evaluate  the  Cauchy  integral  as  described  in  Appendix  B, 
dp/dx  must  be  piecewise  continuous  across  the  midpoints 
x^.  Therefore,  dp/dx  is  taken  as  the  average  of  the  two 
neighboring  grid  point  values  at  every  x^  but  at  those  immedi- 
ately before  and  after  a jump  point,  where  it  is  taken 
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equal  to  the  successive  grid  point  value.  In  this  way 
the  pressure  gradient  discontinuous  jump  is  replaced  by 


an  equivalent  (to  first  order)  linear  variation  across  one 

X 

step  size.  Also,  the  corresponding  value  of  <£. 

XI  X~" 

obtained  by  means  of  equation  (B6)  reproduces  the  same 
singular  behavior  obtained  analytically  for  a discontinuous 
dp/dx  for  Ax  -*•  0 (see  Appendix  E for  details)  . Moreover, 
it  is  regular  for  any  finite  Ax,  thus  allowing  a central 
difference  representation  for  . Numerical  solutions 

obtained  with  such  a model  for  the  pressure  gradient  jump, 
show  first  order  accurate  behavior  versus  Ax  at  the  jump 
points,  thus  confirming  the  validity  of  the  present  approach. 
However,  an  alternate  more  rational  approach  has  also  been 
used  for  further  assessing  the  validity  of  the  present 
numerical  scheme. 


4 . Alternate  Approach  for  the  n+1  Time  Step  Solution 

Most  numerical  schemes  for  solving  the  subsonic  Triple- 
Deck  equations  (Refs.  16  and  17)  use  the  following  inter- 
action law  expression 


7T  r X-£ 


(4.15) 


This  can  be  easily  derived  by  the  one  given  in  equation 

(2.6) 

1 1 P£dC 

D = - i -V-  (2.6) 

xx  tt  * x-£ 

_ nrt  ~ 
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by  simple  integration  by  parts,  which  immediately  gives 


* 


± [£i4L]"  + 1 J-  J P(j)dg  . 

dx  TT  _au  TT  dx  X- 1, 


(4.16) 


At  any  finite  value  of  x the  first  term  in  the  right  hand 
side  of  equation  (4.16)  is  seen  to  identically  vanish  and 
equation  (4.16)  can  be  straightforwardly  integrated  to 
produce  equation  (4.15).  In  order  to  use  equation  (4.15) 
in  the  second  sweep  of  the  present  scheme,  this  is  first 
written  in  finite  difference  form  as 


Dn+1 

l 


-"S3 


Ax 


pd£  _ 
x-C 


CI(xi)  . 


(4.17) 


Note  that  in  equation  (4,17)  the  finite  difference 
representation  for  D^+'*'  is  centered  at  the  midpoint  x^ 
to  avoid  any  differencing  across  the  jump  point.  Also 


note  that  CI(x^)  is  evaluated  without  any  difficulties, 
since  p(£)  is  piecewise  continuous  throughout,  and  only 
once, at  every  location  x^, for  every  n+1  time  level. 
Equation  (4.17)  is  then  differenced  once  to  give 


- Dn+1 
1 + 1 1 

Ax 


Di+1  ' Di-1  CI(x  )-CI(x  ) 

1 )/Ax  = 11  1 


Ax 


Ax 


(4.18) 


which  is  conservative  in  nature  and,  therefore,  exactly 
reproduces  by  integration  the  perfectly  legitimate  equa- 
tion (4.17).  Finally  the  relaxation  integral 
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1-1  „n+l 


x-C 


d£  is  added  to  the  right  hand  side  of 


equation  (4.18)  to  formally  give 


- 2Dn+l  + 

1+1 1 1-1 

Ax2 


CI(xi+1)-CI(x.)  + ^1-1  Dn~t~l_2D*+Dn  ^ 


Ax 


x-? 


(4.19) 


which  is  reduced  to  the  tridiagonal  form  of  equation  (4.5) 
and  inverted  by  Thomas  Algorithm.  The  advantage  of  this 
alternate  approach  is  in  the  evaluation  of  the  Cauchy 
integral  for  the  pressure,  which  has  a regular  behavior 
throughout. 

It  was  found  that  the  length  of  the  region  of  integra- 
tion for  the  pressure  Cauchy  integral  was  extremely  critical 
to  the  convergence  property  of  this  method,  and  if  the 
integration  was  limited  to  the  region  x^,  x^.  as  in  the 
previous  case,  divergent  behavior  was  invariably  observed. 
The  following  analysis  allows  one  to  understand  such 
behavior  and  obtain  a favorable  convergence  property.  In 
the  previous  numerical  scheme  the  interaction  law  was 
effectively  given  by 


1 x1  p£dC 

D = ~ t • 

XX  TT  r X-E, 

X1 

Equation  (4.20)  gives  by  simple  integration  by  parts 

x_ 


(4.20) 


tt:(D„)  = r :rr  (<£  ^ + p(x].)  £n(xI-xi)-p(x1)  £n(xi~x1)  ) 

'C1 


x,  x.-S 
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(4.21) 


Since  p(x^)  = 0 in  the  numerical  scheme,  its  contribution 
to  equation  (4.21)  is  zero.  p(xT)  instead,  even  though 
very  small,  is  ^ 0,  and  neglecting  p(x_j.)  Jinfx^-x^)  in 
equation  (4.21)  is  numerically  equivalent  to  introducing 
a "source  term"  in  the  numerical  procedure,  thus  allowing 
the  solution  to  grow  without  bound,  i.e.,  to  diverge. 

Simply  adding  ^ p(x^)  £n(Xj.-x^)  to  the  value  of  CI(x^) 
given  in  equation  (4.17)  was  found  to  eliminate  the 
problem  and  to  produce  a convergent  behavior  (similar  to 
the  one  typical  of  the  previous  scheme) . This  confirms 
the  foregoing  interpretation  of  the  divergent  behavior 
otherwise  observed.  It  is  believed  that  such  a "source 
term  effect"  is  the  main  cause  of  the  convergence 
difficulties  faced  by  many  schemes  using  the  interaction 
law,  equation  (4.15),  (see  for  example  Ref.  17),  and 
that  a proper  accommodation  of  this  problem  could  signi- 
ficantly improve  their  convergence  property.  For  an 
interesting  example  of  the  importance  of  "source  like  terms" 
to  the  convergence  of  numerical  relaxation  techniques, 
see  Ref.  (31) . 
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V.  RESULTS 


I 


1 . Supersonic  Flow  Results 

The  present  algorithm  as  described  in  Section  III  was 
applied  to  the  case  of  supersonic  flow  past  a parabolic 
hump,  for  which  solutions  of  the  linearized  equations  are 
provided  by  Smith,  (Ref.  8).  The  geometry  of  the  hump  in 
the  nondimens ional  bottom  deck  variables  is  given  in  fig- 
ures (lb  and  2) . The  hump  is  located  in  the  interval 
0 < 0 x < 1,  where  0 is  Lighthill's  constant  (0  = 0.8272) 
(Ref.  8)  . The  hump  profile  is  given  by  hF(x)  = ^x. ( ^ 9 x ) 

U 

in  the  interval  0 £ 0 x <_  1;  its  maximum  height  is  — =• 

40 

corresponding  to  0x  = 0.5.  A small  hump  height  (h  = 0.1) 
was  first  taken  into  consideration  in  order  to  check  the 
present  scheme  by  direct  comparison  with  Smith's  linear 
results  (Ref.  8).  Figure  (2)  provides  a comparison  between 
the  normalized  wall  shear  and  pressure  obtained  by  the 
present  numerical  scheme  for  the  nonlinear  set  of  equations 
and  Smith's  linear  solutions.  The  agreement  is  quite 
convincing  and  no  other  comparisons  were  sought.  Step  size 
studies  have  been  performed  in  order  to  further  confirm  the 
reliability  of  the  present  method  as  well  as  verify  its 
accuracy.  Figures  (3)  and  (4)  present  step  size  studies 
for  normalized  pressure  and  wall  shear  at  the  two  most 
critical  points,  that  is  the  leading  and  trailing  edges  of 
the  hump  (0x  = 0 and  0x  = 1) . As  the  results  clearly  in- 
dicate, the  scheme  is  verified  to  be  first  order  accurate 
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in  the  longitudinal  x direction  and  second  order  accurate 
in  the  normal  z direction.  After  that  the  reliability  of 
the  scheme  was  satisfactorily  assessed  by  the  comparisons 
and  step  size  studies  just  presented,  attention  was  devoted 
to  hump  heights  of  order  one.  In  this  case,  since  the  full 
nonlinear  set  of  equations  only  holds,  a numerical  solution 
is  the  only  feasible  one  and  is  also  the  only  means  for 
assessing  the  occurrence  of  separation.  Numerical  results 
were  obtained  for  values  of  h up  to  6,  for  which  separation 
and  reattachment  occur  around  the  trailing  edge  of  the  hump. 
Figures  (5),  (6)  and  (7)  present  normalized  wall  shear, 

pressure  and  displacement  thickness  for  values  of  h up  to 
4,  compared  with  those  corresponding  to  h = 0.1,  in  order 
to  assess  the  influence  of  the  nonlinearities  on  the  solu- 
tion profiles.  As  appears  in  figures  (5)  and  (6)  the 
effect  of  nonlinearity  is  seen  to  be  as  important  for  h 
increasing  from  2 to  4 as  for  h increasing  from  0.1  to  2, 
(which  is  reasonable)  and  the  maximum  deviations  from  the 
linear  behavior  occur  in  correspondence  to  the  hump  itself. 
Figure  (7)  shows  the  strange  behavior  (already  anticipated 
in  Section  III)  taken  by  the  displacement  function,  which, 
as  the  total  displacement  (6  = D for  0x  > 1),  tends  to  a 
nonzero  asymptotic  value  as  x -*•  °°.  Such  behavior  is  due 
to  truncation  error  and  is  a considerable  entity  since  the 
scheme  is  only  first  order  accurate  in  Ax.  The  insert  of 
Figure  (7)  shows  that  at  0x  = 2 the  value  of  6 extrapolated 
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to  zero  step  size  Ax  is  positive  (not  negative  as  appears 
in  figure  (7)  and  is  therefore  decaying  as  x -*■  00  as  it 
should.  Fortunately  this  improper  behavior  of  5 for  finite 
fairly  large  values  of  Ax  does  not  influence  the  pressure 
and  shear  distributions  which  are  related  to  the  derivatives 
of  the  total  displacement  and  are  therefore  not  affected  by 
the  constant  downstream  shift  of  6.  The  impact  of  non- 
linearity is  also  shown  for  the  velocity  profiles  at  three 
locations,  leading  edge,  maximum  height  and  trailing  edge 
of  the  hump  (9x  = 0,  0.5  and  1).  The  (linear  theory 

(Ref.  8)  normalized  longitudinal  velocity)  profiles  are 
presented  here  for  convenience  in  figure  (8)  for  two  values 
of  h.  The  flow  field  is  slightly  decelerated  at  9x  = 0, 
strongly  accelerated  at  9x  = 0.5  and  strongly  decelerated 
near  the  wall  at  9x  =1,  where,  for  h = 4,  the  wall  shear 
appears  to  be  already  negative.  This  is  in  agreement  with 
physical  expectation  as  well  as  linearized  equations  results. 

Finally  the  total  shear  is  shown  in  figure  (9)  to  better 
display  the  large  gradients  and  the  wide  variations  for  the 
shear  which  is  seen  to  become  negative,  already  for  h = 4 
and  more  so  for  h = 6, around  the  trailing  edge  of  the  hump. 
This  confirms  Smith's  assumption,  based  on  linear  results, 
that  the  most  likely  place  for  separation  to  occur  would 
indeed  be  at  the  backside  of  the  hump  (Ref.  8) . Figure 
(9)  also  confirms  the  capability  of  the  present  scheme  to 
properly  assess  separated  flows.  The  results  shown  in 
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figures  (4)  through  (9)  have  been  obtained  with  a rather 
course  mesh  in  both  the  x and  z directions.  They  are  be- 
lieved to  be  qualitatively  correct  but  only  step  size  studies 
could  reasonably  assess  their  accuracy. 

In  the  present  study  the  outer  edge  boundary  conditions 
were  imposed  at  z = 10  for  the  smaller  h values  and  z = 15 
for  the  higher  ones.  In  all  cases  it  was  numerically  veri- 
fied that  the  u velocity  profile  became  linear  (u-*z-6) 
over  several  mesh  points  below  the  outer  edge.  The  upstream 
initial  and  downstream  boundary  conditions  were  set  respect- 
ively at  x = - 5 and  x = 15  and  were  numerically  verified 
to  be  at  a sufficient  distance,  since  further  increase  of 
the  region  of  integration  produced  no  appreciable  change  in 
the  results.  The  numerical  procedure  was  found  to  coverge 
(average  absolute  error  on  the  total  displacement  height, 

[£  |D^+1  - D? | ]/(I-l),  < 10  within  20  iterations  (for 
h up  to  1)  and  the  computer  time  required  for  a typical 
case  (0Ax  = 0.1  Az  = 0 . 25,  corresponding  to  201  by  41  element 
matrices  for  the  u and  w velocities)  was  about  15  CPU  seconds 
(AMDAHL  470  double  precision  arithmetic) . The  approximate 
optimal  time  step  for  convergence  was  numerically  found  to 
be  given  by  A = 1 . In  the  presence  of  separation,  convergence 
became  slower  and  X had  to  be  increased  up  to  6 (for  h = 6) 
to  avoid  divergent  behavior.  For  the  present  geometry  it 
is  felt  that  h = 6 is  close  to  the  highest  value  of  h for 
which  the  present  scheme  could  converge  because  of  the 
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horrendous  gradients  due  to  the  geometry  itself.  The  same 
scheme  and  a totally  second  order  version  of  the  same  were 
also  found  to  be  very  efficient  for  supersonic  flow  past  a 
wedge.  In  this  case  convergence  was  easily  obtained  for 
normalized  wedge  angles  up  to  3,  corresponding  to  a fairly 
large  separation  bubble  and  the  results  agreed  very  well 
with  those  of  Refs.  (13-15). 

At  the  intial  time  level  the  pressure  gradient  and 
displacement  thickness  were  taken  equal  to  zero  at  all  x^ 
locations  in  order  to  start  the  relaxation  procedure. 

This  physically  corresponds  to  the  introduction  of  the  dis- 
turbing hump  as  time  equals  zero  in  an  otherwise  undisturbed 
uniform  shear  flow.  For  the  higher  hump  cases,  h was  in- 
creased in  gradual  steps  at  every  new  time  level  tn  until 
the  final  value  was  reached.  This  was  done  in  order  to 
avoid  divergent  behavior  produced  by  an  otherwise  too  large 
time  rate  of  change  in  the  transient  solution,  which  as 
already  mentioned  has  no  physical  meaning  anyway. 

2 . Subsonic  Flow  Results 

The  subsonic  version  of  the  present  scheme  as  described 
in  Section  IV  was  applied  to  the  solution  of  flow  past  the 
parabolic  hump  of  figures  (lb  and  2)  as  well  as  the  quartic 
hump  of  figure  (lc) . For  the  first  case  linearized  equation 
results  are  provided  by  Smith  (Ref.  8),  whereas  a numerical 
solution  corresponding  to  a fairly  large  separation  bubble 
is  given  by  Sykes  (Ref.  23)  in  the  second  case.  The  quartic 
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hump  of  figure  (lc),  whose  profile  is  given  by  hF(x)  = 

2 ^ 

h(l-x  ) for  - 1 <_  x <_  1,  has  also  the  advantage  of  having 
continuous  slope  at  its  leading  (x  = - 1)  and  trailing 
(x  = 1)  edges  and  therefore  does  not  require  any  special 
numerical  accommodation  of  the  present  algorithm  at  those 
locations.  In  both  cases  numerical  solutions  have  been 
obtained  by  means  of  either  numerical  scheme  discussed  in 
Section  IV. 

Subsonic  flow  past  the  parabolic  hump  of  figure  (2) 
will  be  considered  first.  The  small  hump  height,  h = 0.1, 
has  been  again  considered  as  a test  case  in  order  to  compare 
the  present  numerical  results  with  Smith's  linear  ones  (Ref. 
8) . Figures  (10)  and  (11)  show  the  normalized  wall  shear 
and  pressure  profiles  presently  obtained  compared  with  those 
of  Ref.  (8).  Whereas  a reasonable  agreement  is  obtained  for 
the  case  of  the  shear,  the  two  pressure  distributions  show 
disturbing  differences  especially  immediately  upstream  and 
downstream  of  the  hump.  Smith  has  kindly  recalculated  the 
pressure  and  shear  profiles  and  his  "improved"  results  are 
now  indistinguishable  from  the  present  ones  (Ref.  30). 

The  present  results  given  in  figures  (10)  and  (11)  have 
been  obtained  by  means  of  the  second  version  of  the  algorithm 
(scheme  B)  described  in  Section  IV.  A longitudinal  step 
size  study  has  been  performed  in  order  to  assess  the  accuracy 
and  reliability  of  the  solution  technique.  The  results  are 
given  in  figures  (12)  and  (13)  where  wall  shear  and  pressure 
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are  seen  to  be  first  order  accurate  versus  Ax  at  three 
different  locations  (0x  = 0,  0.5,  1).  Figure  (13)  is  the 
more  interesting  in  that  it  shows  2 different  pressure 
results  for  each  location  x.  The  circles  indicate  the 
pressure  values  obtained  by  simple  Euler  quadrature  of  the 
pressure  gradient  distribution  at  the  final  time  step  n+1, 
after  convergence  has  been  reached.  The  squares  indicate 
the  values  obtained  by  a second  approach,  exploiting  the 
reciprocity  property  of  the  Riemann  integrals.  In  fact, 
besides  the  interaction  law  equation  (4.15) 


1 , P£dC 

Dx  - ? * 


(4.15) 


the  following  conjugate  integral  expression  is  also  true, 


p = - = & d£ 

* ir  * x-£  * 


(5.1) 


Therefore,  after  the  converged  solution  has  been  reach- 
ed by  means  of  the  present  numerical  scheme,  the  pressure 
was  evaluated  by  means  of  equation  (5.1).  The  Cauchy  in- 
tegral (5.1)  was  evaluated  by  means  of  equation  (B6)  where 
Dj.  was  replaced  by  a backward  finite  difference  represent- 
ation and  the  correction  was  given  by  a central  finite 

difference  representation.  The  two  different  pressure 
evaluations  are  shown  in  figure  (13)  to  be  both  first  order 
accurate  and  at  locations  0x  = 0 as  well  as  0x  = 0.5  they 
tend  to  the  same  limit  values  as  Ax  goes  to  zero.  For 
0x  = 1 the  extrapolated  results  are  slightly  different  and 
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this  is  probably  due  to  the  finiteness  of  the  integration 
region  for  the  Cauchy  integrals  as  well  as  to  the  growth 
of  the  Goldstein  layer  originating  at  the  leading  edge  of 
the  hump  (Ref.  15) . 

The  effect  of  the  finite  Az  step  was  also  investigated 
but  it  was  found  to  be  negligible  at  the  level  of  Az  = 0.2 
used  in  all  present  calculations.  (It  is  worth  remembering 
that  the  scheme  is  second  order  accurate  in  the  normal  z 
direction.) 

The  outer  edge  boundary  conditions  were  imposed  at 
z = 8 to  15  with  increasing  h.  The  range  of  integration 
in  the  longitudinal  direction  was  taken  to  be  -5  <_  0x  < 8 
and  moving  the  range  further  upstream  and  downstream  the 
two  limits  was  seen  to  produce  minor  changes  in  the  solution. 
Also,  introducing  the  downstream  far  field  asymptotic  eval- 
uation of  the  Cauchy  integral  residual  as  described  in 
Appendix  D was  seen  to  modify  the  results  by  about  two  per- 
cent or  less.  Optimizing  the  procedure  could  probably  have 
produced  the  same  level  of  accuracy  with  a shorter  numerical 
integration  region  but  this  aspect  was  not  pursued  further 
at  the  present  time.  All  the  results  presented  here  have 
the  Cauchy  integral  region  of  integration  limited  to  the 
range  - 5 0x  < 8 . 

The  approximate  optimal  time  step  for  the  subsonic 
algorithm  was  obtained  for  X = 2 but  again  it  had  to  be 
increased  for  the  high  h cases  in  order  to  avoid  divergent 
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behavior.  All  the  results  presented  here  have  been  obtained 
by  imposing  a limit  of  ten  to  the  number  of  iterations  in 
the  second  sweep.  The  total  procedure  converged  (average 
error  of  the  total  displacement  < 10  ) in  about  50  to  80 

time  cycles  for  h £ 2,  depending  on  the  value  of  the  step 
size  Ax  (for  lower  values  of  Ax  more  iterations  were  needed) . 

A typical  calculation  required  about  two  to  three  minutes 
of  computer  time  (AMDAHL  470,  double  precision  arithmetic) 
a very  small  amount  for  a subsonic  Triple-Deck  solver. 
Convergence  became  very  difficult  for  h > 4 and  the  need  of 
further  assessing  the  capability  of  the  present  scheme  of 
handling  reverse  flow  conditions  irrespective  of  slope  dis- 
continuity singularities  led  to  consideration  of  Sykes 
quartic  hump  (Ref.  23) . The  smooth  nature  of  such  a geometry 
(see  figure  lc)  also  allowed  the  use  of  both  schemes  described 
in  Section  IV  (and  called  schemes  A and  B for  convenience) 
to  provide  a double  check  of  the  results.  This  was  again 
obtained  by  means  of  a longitudinal  step  size  study.  Fig- 
ure (16)  shows  the  wall  shear  results  obtained  by  both 
schemes  at  the  3 locations  x = - 1,  0,  1,  for  various  values 
of  Ax.  The  extrapolated  results  basically  coincide  to 
provide  further  credibility  to  the  present  approach.  Figures 
(17)  and  (18)  show  the  pressure  values  obtained  by  Euler's 
quadrature  of  the  pressure  gradient  and  by  evaluation  of 
the  inverse  Cauchy  integral,  equation  (5.1),  using  schemes 
A and  B respectively.  Again  all  the  extrapolated  results 
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agree  very  well  to  prove  once  more  the  reliability  of  the 
numerical  procedure.  Figure  (19)  provides  normalized  wall 
shear  and  pressure  distributions  obtained  for  h = 0.1  by 
means  of  the  first  version  of  the  algorithm,  scheme  A. 

For  this  case  the  outer  boundary  conditions  were  imposed 
at  z = 10  and  the  region  of  integration  was  taken  to  be 
- 6 £ x <_  10 . 

Finally  the  h = 3 case  was  considered  in  order  to 
test  the  ability  of  the  present  scheme  to  handle  the  con- 
siderable reverse  flow  bubble  anticipated  by  Sykes  numer- 

: 

ical  predictions  (Ref.  23) . For  this  case  convergence  was 
very  slow  even  when  the  Reyhner  Fliigge-Lotz  approximation 
(Ref.  32)  was  used  and  the  same  phenomenon  observed  in 
Ref.  (11)  for  the  highest  Reynolds  number  cases  was  en- 
countered. The  solution,  while  steadily  converging  every- 
where else,  presented  an  oscillatory  behavior  in  time  near 
the  reattachment  point,  probably  due  to  a source  effect 
induced  by  the  switching  of  the  finite  difference  repre- 
sentation (forward-backward)  of  the  longitudinal  convective 
term.  Upwinding  the  normal  convective  term  was  found  some- 
what helpful,  and  a way  for  completely  eliminating  such  an 
oscillatory  behavior  around  reattachment,  was  to  set  both 

convective  terms  equal  to  zero  wherever  the  coefficient  u 

* *o 

(here  taken  as  (u^+U2  )/2)  in  the  longitudinal  convective 
terms  uux  of  momentum  equation  (3. 10) was  found  to  be  negative. 

Figures  (20)  and  (21)  show  the  shear  and  pressure  profiles 
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obtained  by  introducing  this  approximation  into  the  present 
scheme,  as  compared  to  Sykes  numerical  results  (Ref.  23). 

The  agreement  is  satisfactory  if  we  c6nsider  that  Ax  = 0.2 
in  the  present  first  order  accurate  scheme  whereas,  Ax  = 1/12 
in  Sykes  second  order  accurate  scheme . More  accurate  solu- 
tions could  not  be  obtained  in  this  study  because  divergent 
behavior  was  observed  for  smaller  values  of  the  step  size 
Ax.  Such  phenomenon  is  exactly  equivalent  to  what  was 
already  encountered  in  high  Re  separated  flow  solutions  of 
the  interacting  boundary  layer  equations  in  the  vorticity 
stream  function  formulation  (Ref.  11) . Methods  for  correct- 
ing this  divergent  behavior  are  currently  under  investigation. 

Once  the  source  of  such  "instabilities" is  fully  understood 
it  is  believed  that  accurate  solutions  for  "high  Reynolds 
number"  separated  flows  will  be  feasible  by  means  of  the 
present  Triple-Deck  solver.  Understanding  this  effect  is 
also  important  in  numerical  solution  of  the  Navier  Stokes 
or  interacting  boundary  layer  equations  since  we  have  en- 
countered the  same  difficulty  there  as  well.  • ' 
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APPENDIX  A 


INVERSION  PROCEUDRE  FOR  THE  * SWEEP  FINITE 


DIFFERENCE  EQUATIONS 

The  two  sets  of  infinite  difference  equations  given  by 
equations  (3.11),  (3.12)  and  (3.13),  (3.14)  are  easily  re- 
cognized to  be ’of  the  form 


w.  =w.  . -c.u.  - c . u.  . +d. 
3 3*1  33  3 3-1  3 


Ajuj-i  + Vj  + c:Vi  * Hj  ' a5j“j  (A2: 

where  the  coefficients  c^  through  a5^  can  be  immediately 
identified  by  direct  comparison  with  equations  (3.11) 
through  (3.14).  Such  a coupled  set  of  finite  difference 
equations  can  be  directly  inverted  (see  Ref.  22)  by  means 
of  the  simple  recursion  relation 

uj " Vj-i + f j + Vj-i  <A31 

where  the  coefficients  E^,  F ^ are  given  at  location  j 
by  the  following  equations. 


- Aj  - <cjGi+i+a5i)a- 


[V  cjEj+r(cjGj+i+a5j,ej1 


H.  - C .F  . . (C  . G . . +a5  . ) d . 
3 3 3+1  3 3+1  3 3 


(Bj  + C.Ej+1(C.Gj+1+a5.)a.] 


The  complete  solution  for  the  two  arrays  and  can  be 
directly  obtained  by  means  of  equations  (Al)  and  (A3)  once 
the  wall  values  u^  and  w^  are  known  together  with  all  the 


coefficients  E.,  F.,  G..  The  values  u.  and  w.  are  directly 
3 3 3 1 1 

provided  by  the  wall  boundary  conditions  which  in  this  case 
are  u^  = w^  = 0 (see  equations  (3.16)  for  comparison).  The 
coefficients  E ^ , FV,  G ^ can  be  obtained  by  means  of  equations 
(A4)  through  (A6)  once  their  values  at  the  last  grid  point 
Z_  are  known.  These  are  determined  by  means  of  the  outer 
edge  boundary  condition  equations  (3.18)  or  (3.19)  as  des- 
cribed below.  The  last  three  points  in  the  computational 
mesh  are  indicated  for  simplicity  as  3,  2,  1 (see  figure 
Al)  and  the  following  second  order  representations  of  the 
longitudinal  velocity  u at  points  2 and  1 are  written  by 
means  of  Taylor  series  expansions  : 


Figure  Al . Boundary  Condition  Treatment. 


60 


' Az2  ' 


and 


u2  = u3  - Azu3  + ~y~  u3  + 0 ( Az  ) 


' 2 " -5 

= u3  - 2Azu3  +2Az  u3  + 0(Az  ) 


(A7 ) 


(A3) 


(Note  that  ’ indicates  here  partial  differentiation  with 

I 

respect  to  z.)  By  replacing  u3  by  its  known  value  x 

It 

(x  = 0 or  1)  and  eliminating  u3  in  equations  (A7)  and  (A8), 
the  following  second  order  accurate  relation  for  the  values 
of  u at  the  three  locations  of  figure  Al  is  obtained 


u^  = 4u2  - 3u3  + 2xAz 

The  finite  difference  equation  (A2)  is  then  written  at 
grid  point  2 to  give 


(A9 ) 


A2u^  + B2u2  + *-2u3  = H2  “ a^2w2 


(A10) 

which,  after  eliminating  u^  by  means  of  equation  (A9),  can 
be  written  as 


B„-4A„  H2~2xA2Az  a52 


2 2 

U3  ' C2-3A2  U2  + C2-3A2 


C2-3A2 


(All) 


Note  that  equation  (All)  is  nothing  else  but  the  recursion 
formula  (A3) , written  at  the  last  grid  point,  3 = J,  and 
provides  the  sought  after  relations  for  the  coefficients 
E3 , F3,  G3 , namely 


E ..W 

j C_-3A_ 


(A12 ) 
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Appendix  is  performed  by  a subroutine  which  receives  from 


the  main  program  all  the  arrays  c^  through  a5^  and  the 
appropriate  value  of  X (0  or  1)  and  returns  to  the  main 


r 


APPENDIX  B 


NUMERICAL  SOLUTION  OF  CAUCHY  INTEGRALS 


The  problem  considered  here  is  the  numerical  evaluation 
of  a Cauchy  integral  of  the  type 
b 


CI(x)  = f 


f (€)dg 
x-C 


for  a < x < b 


(Bl) 


where  f(£)  is  a piecewise  continuous  function  in  the 
interval  a,b. 

In  order  to  numerically  evaluate  CI(x),  the  interval 
of  integration  a,b  is  divided  into  an  arbitrary  number  of 
equally  spaced  grid  points  x^  (i=l,  2,  ...  I)  such  that 

x^  = a and  x^  = b.  The  midpoints  of  each  grid  are  also 

• • • — — 
important  and  are  indicated  as  x^  (x^  = x^  + , 

i ■ 1,  2,  ...  1-1) . 

CI(x)  can  be  numerically  evaluated  at  any  location  x^ 
by  the  following  expression. 


!~1  _ Xj+1  1-1  i x.-x. 

CI(x.)  = l f ( x . ) f ~-  = l f(x.)  y ^n(ri — 1 
j = l 11  Xj  x^S  j=l  : xi"xj 


-) 


j+1 

(B2) 


or,  to  the  same  order  of  accuracy, 

1-1  f(x.)+f(x.  .)  . X.-X.  2 

CI(x,l  - l J-5 ItL  1 tn  (_x—  3-) 

3=1 


(B3) 


Note  that  equations  (B2)  or  (B3)  correspond  to  the  simple 
trapezoidal  rule  in  which  the  function  f(£)  has  been 
averaged  over  each  mesh  of  the  integration  domain  and  the 
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remaining  term  d£/(x-£)  has  been  analytically  integrated 
over  each  grid  to  properly  account  for  the  Cauchy  principal 
value  nature  of  the  integral  CI(x).  Therefore,  the 
accuracy  of  equations  (B2)  and  ( B 3 ) might  be  expected  to 
be  second  order  versus  Ax.  However,  application  of 
equation  ( B 3 ) to  the  following  model  problem, 

1 -2  . 

CI(x)  = £ dC  = x(x<ln  " 2)  (B4) 

(Ref.  25),  produced  only  first  order  accurate  solutions. 

Figure  22  clearly  shows  that  Cl (0.45)  as  well  as  Cl (-0.15) 

evaluated  by  means  of  equation  (B3)  linearly  approach  the 

exact  analytical  values  (in  the  limit  of  zero  step  size.  Ax). 

In  general,  a first  order  accurate  numerical  integration 

is  able  to  exactly  reproduce  the  analytical  value  of  the 

sought  after  integral  only  for  a constant  value  of  the 

integrand  (Ref.  24) . In  the  present  special  case  then, 

equation  (B2)  should  be  able  to  provide  the  exact  solution 

for  ^5'"  r*  only  if  f{£)  is  a constant  in  any  interval 
a x~^ 

x.-x . . . This  was  numerically  verified  for  the  following 
J J 

model  problem,  £ = xHn  - 2,  (Ref.  25)  . 

— 2 -L“X 

The  numerical  integration  scheme  provided  by  equation  (B2) 
or  (B3)  as  well,  is  therefore  only  first  order  accurate 
and  such  an  apparently  anomalous  result  needs  further 
analysis.  A first  understanding  of  such  behavior  can  be 
obtained  by  noting  that  in  equation  (B2)  or  (B3)  the 
integration  across  the  step  x^-xi+^  (in  which  the  argument 
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of  CI(x^)  becomes  singular)  gives  zero  contribution  to  the 
integral  itself.  However,  the  exact  contribution  to  the 
total  integral  CI(x^)  due  to  the  integration  across  such 
an  interval  is  zero  indeed  only  if  the  function  f(£)  is 
even  around  x^,  the  midpoint  of  the  interval  itself. 

In  general,  instead,  it  might  be  reasonably  expected  to 
be  of  order  Ax,  the  length  of  the  interval  itself. 

Neglecting  such  contribution  could  therefore  be  sufficient 
to  produce  the  first  order  accurate  behavior  of  the  present 
scheme.  The  above  argument  received  numerical  support 
by  the  solution  of  the  following  integral 

ci,(*)  = # i , IB5) 

j.  _ ^ 

• i 

obviously  chosen  because,  in  light  of  the  foregoing  dis- 
cussion,  equation  (B2)  or  (B3)  should  provide  a second 
order  accurate  solution  for  CI^(0.45)  whereas  only  a 
first  order  accurate  one  for  CI^Cx)  at  any  other  location 
x 7*  0.45.  Figure  23  shows  that  this  is  indeed  the  case, 
thus  supporting  the  previous  discussion.  A complete 
understanding  of  the  problem  will  be  shortly  obtained 
by  analyzing  a general  second  order  accurate  integration 
scheme  for  equation  (Bl) , provided  below.  Such  a second 
order  scheme  should  be  able  to  resolve  the  previously 
discussed  anomaly  and  also  to  reproduce  the  exact  analytical 
value  of  CI(x)  for  any  piecewise  linearly  varying  function 
f(5).  This  can  be  easily  accomplished  by  extracting  the 
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average  values  of  the  function  f(£)  and  its  first  derivative 
over  every  integration  mesh  interval,  analytically 
integrating  the  remaining  part  of  the  integrand  to  properly 
account  for  its  singular  nature,  and  finally  summing  all 
these  contributions  in  the  following  way: 
b 


f <5)d€ 


f(M  1-1  j+1 

CI(xi)  = t H 

a x^S  j=l  Xj  x^S 


1-1  f(x.  _)+f(x.)  Xj+1  1-1  f(x.  ,)-f(x.)  Xj+1  £-x. 

- I — ^ l -1 £ L< 


j-1 


xj  V5  i’1 


x.  xL-( 


- i-i 

3-1  xi-xj+i 


y ri.  ( *3+1  dJL 

A Ai 


j = l 


Ax 


^ f'vV  £ |=f  - / d^; 


Xj 


Xj 


finally. 


1-1  f (x,xl ) +f (x.. ) 

l 

j-1 


x.-x.  2 


CI(xi)  = l ^2 J Zn  — 3 > 


Xi~Xj+l 


1-1  f(x.  ,)-f(x.)  , X.-X.  2 

l 3-  t(x  -x  ) J ) - Ax]  . 

*i  si  J y »y 

3 i j+1  (B6) 


Note  that  equation  (B6)  essentially  contains  the  same 
term  as  equation  (B3)  plus  corrections  proportional  to 


the  local  slope  of  the  function  f(5),  [f (x j+i) -f (Xj ) ] /Ax  = 


f (Xj)+0(Ax  ).  Such  a contribution,  for  the  case  (i-j)  / 1, 


can  be  written  as 


66 


f'  (Xj)  [(i-j)Ax  in  2(i~rjf+f  ■ Ax]  + 0(Ax2)  (b7) 

and  is  easily  seen  to  be  a first  order  correction.  Equations 
(B6)  and  (B7)  also  show  that  if  f(£)  is  even  ( f ' (5)  odd) 
around  the  x^  location,  the  first  order  terms  exactly 
cancel  each  other,  thus  producing  the  second  order  accurate 
behavior  observed  in  Figure  23.  Equation  (B7)  further 
shows  that  the  dominant  0(Ax)  contributor  is  due  to  the 
step  size  containing  the  singularity, as  previously 
anticipated.  In  fact,  for  (i-j)  ? 1,  Equation  (B7)  gives, 


by  Taylor  series  expansion  in  terms  of  the  "small 
parameter"  l/(i-j), 

f'(Xj)  [Ax(j^4-  + H.O.T. ) - Ax]  + 0 (Ax2)  , (B8) 

which  is  seen  to  be  a "small"  0(Ax)  contribution. 

For  i = j.  Equation  (B7)  gives  instead, 

(-f'fXj)  + H.O.T.)  Ax  + 0 ( Ax2 ) , (B9 ) 

which  is  the  main  first  order  contribution  neglected  by 
equation  (B2)  or  (B3) , whose  first  order  accurate  behavior 
is  now  completely  understood. 


It  is  important  to  mention  that  in  equation  (B6)  the 


group 


f (x.)+f  (x.+1) 


can  be  replaced  by  f(x^)  and  the  one 


f(x  )-f(x.) 


f(x  )-f(x.  x) 


, to  produce  four  second 


order  accurate  schemes.  Here  two  of  them,  namely  the  one 
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given  in  equation  (B6)  and  that  obtained  with  both  the 

above  indicated  modifications,  have  been  applied  to  the 

model  problem,  equation  (3),  and  are  shown  to  recover 

the  exact  solution  with  second  order  accurate  behavior 

in  figure  24.  They  also  have  been  able  to  reproduce 

the  exact  analytical  solution  for  the  problem, 

1 . 

<£  d£  = x£n  -z — — - 2,  (Ref.  25)  , for  any  of  the  step 
h,  l— x 

sizes  used  for  the  results  of  figure  24,  thus  confirming 
once  more  their  second  order  accurate  nature. 

It  must  be  acknowledged  that  the  present  technique 
is  a special  application  of  more  general  methods  for 
computing  two  dimensional  potential  flows  about  arbitrary 
bodies  by  means  of  source  and  or  vortex  singularity 
distribution  (see  Refs.  26,  27,  28).  Dilley,  for  example, 
has  already  used  generalized  expressions  of  equation  (B6) 
(Ref.  29) . His  failure  of  achieving  second  order  accuracy 
is  most  probably  due  to  the  curved  nature  of  the  bodies 
considered  in  his  study  or  to  the  variable  step  size  of 
his  computational  domain. 
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APPENDIX  C 


A VECTORIZATION  OF  LOGARITHMS  IN  THE 


CAUCHY  INTEGRAL  SOLVER 


It  has  been  seen  in  Appendix  B that  an  extremely 
large  number  of  logarithms  must  be  computed  in  order  to 
evaluate  the  Cauchy  integrals  of  the  subsonic  interaction 
law,  equation  (2.6).  Such  logarithms  are  of  the  form 

AL.  . = j Hn  (_-~  j ) (Cl) 

Xi'Xj+l 


with  2 <_  i <_  1-1,  2 <_  j <_  1-1,  and  could  be  easily  cal- 
culated once  and  for  all  and  stored  in  a square  matrix 
of  order  1-2,  to  avoid  the  extremely  large  computational 
effort  of  evaluating  them,  when  they  are  needed.  It  is 
worth  mentioning  that  for  every  iteration  in  the  second 
sweep  of  the  numerical  procedure  described  in  Section  IV, 

the  single  operation  of  equation  (Cl)  is  performed 
2 

(1-2)  times.  Considering  that  such  second  sweep  cal- 
culations can  be  performed  up  to  a 1000  times  (10  * 100) 
and  that  logarithmic  calculations  are  extremely  time 
consuming  on  the  computer,  it  is  easy  to  understand  the 
great  need  of  calculating  equation  (Cl)  once  and  for  all. 
The  use  of  a square  matrix  though  requires  large  storage 
for  small  step  size  computations.  It  is  sufficient  to 
mention  that  for  (1-2)  = 200,  a double  precision  evaluation 


of  AL^  j would  require  320  Kilobytes  of  central  memory, 
which  would  rise  up  to  the  prohibitive  figure  of  2000 
for  (1-2)  = 500.  Fortunately,  a simple  vector  of  21-3 
length  can  serve  the  purpose  equally  well,  if  a constant 
grid  is  used  in  the  computational  scheme.  It  is 
sufficient  to  observe  that  the  argument  of  equation  (Cl) 
only  depends  on  the  relative  distance  between  the 
points  x^  and  x^  (see  figure  Cl) . 
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j > i 
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■*  X >-■  X > X i X- 
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x . 
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Figure  Cl.  Grid  Distribution. 


If  (i-j)  < 0 


AL.  . = AL. 
1,3  k 


r (j-i)  Ax-Ax/2, 
*n  (j-i)  Ax+Ax/2 


£n  2 ( j - i) -1  (C2) 

2 (j-i)  +1  ' 


if  (i-j)  > 0 


ALi,j  - “k  - tn  friffir  ■ 


( C 3 ) 


and  finally,  if  (i-j)  = 0 


AL.  . = AL,  = 0 
1,3  * 


(C4 ) 


It  is  sufficient  to  calculate  once  and  for  all  the 
values  of  equations  (C2)  through  (C4)  and  store  them  in 
the  vector  AL^.  For  k = I,  AL  = 0 and  for  1+1  < k < 21-1, 
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1 


AL^  will  contain  the  values  corresponding  to  (i-j)  > 0, 

whereas  for  1 < k < I,  AL^  will  contain  the  values 

corresponding  to  (i-j)  < 0.  It  is  sufficient  to  define 

k = I + i-j,  and  the  correct  value  (AL  =AL.  .)  will  be 

K 1 / 3 

immediately  available  whenever  necessary.  Fortran 
programming  of  this  procedure  is  straightforward  and 
very  short.  Its  use  has  produced  considerable  saving  of 
computing  effort.  (For  a case  in  which  2 minutes  of 
AMDAHL  470  are  presently  used,  repetitive  calculation  of 

II 

the  logarithms  led  to  a total  computational  effort  of 
10  minutes.) 

Finally,  it  must  be  noted  that  the  vector  AL^  is 
symmetric  around  its  central  element  AL^  and  properly 
exploiting  such  a property  could  lead  to  a further 
reduction  of  the  computer  storage  from  21-3  to  just  1-1 
double  precision  locations.  Considering  though  the  more 
elaborate  logic  and  the  larger  computational  effort 
needed  to  implement  such  a procedure,  the  present  approach 
is  definitely  the  simplest  and  most  efficient  one. 
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APPENDIX  D 


A 


ANALYTICAL  EVALUATION  OF  FAR  FIELD  ASYMPTOTIC 
CONTRIBUTION  TO  THE  CAUCHY  INTEGRALS 

The  problem  of  interest  here  is  the  evaluation  of  an 
integral  of  the  type 

X1 

/ for  > (Dl) 

7 for  x < x . ( D2 ) 

XI 

Note  that  in  equations  (Dl)  and  (D2) , the  Cauchy  principal 
value  sign  is  missing  because  the  argument  becomes  singular 
outside  the  interval  of  integration.  An  accurate  evaluation 
of  the  integrals  (Dl)  and  (D2)  is  desirable,  due  to  the 
slow  algebraic  decay  of  f (£)  as  ( + +«  or  as  £ -»■  -00 
in  most  Triple-Deck  equation  applications.  In  this  appendix 
only  the  solution  of  equation  (D2)  will  be  considered, 
for  a function  f(£)  decaying,  as  £ •+■  + °°,  with  a negative 
integer  power  law.  This,  because  for  subsonic  flow  past 
a hump  the  pressure  gradient  asymptotic  decay  is  given 
(Ref.  8)  as: 

§£  = f m = k C2  (for  C » 1)  . (D3) 

Extension  to  larger  values  of  the  exponent  is  straight- 
forward. The  particular  integral  of  interest  is  then 
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00 


(D4) 


I(x) 


k / 

x 


dC 

(x-£K3 


A simple  change  of  variable  of  integration  produces 


I(x) 


^ / 

X X.  /X 


dn 

(1-n) n3 


(D5) 


and  the  integral  of  equation  (D5)  can  be  solved  by  the 


method  of  partial  fractions.  In  fact,  since 

1 1,11.1 


n (1-n) 


l_  + + I + 

n3  n2  n i-' 


I(X)  = Kr  [——*  - i + s.n  n - 4n(n-l)]  • (D6) 

xJ  2n  n x /x 

Equation  (D6)  is  only  valid  for  x > 0,  and  a very  similar 
one  should  be  derived  for  x < 0.  Simple  regrouping  of 
the  logarithmic  terms,  though,  produces  the  following  final 
expression  for  I (x),  which  is  easily  verified  to  be  valid 
for  negative  as  well  as  positive  values  of  x 


I(x) 

k 


[In  ( 


(*-)  1 
XI 


(D7) 


Note  that  Xj  is  the  lower  limit  of  integration  and 
has  to  be  large  enough  to  guarantee  a satisfactory  accuracy 
for  the  asymptotic  representation  of  f(£),  equation  (D3) . 
x,  instead,  is  the  variable  point  at  which  the  integral 
r (x)  is  evaluated  and  can  vary  within  the  interval 
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x^  < x < x . A special  treatment  is  of  course  required 
for  x = 0.  In  this  case  the  integration  trivially 
produces 

GO 

I(x)/k  = / = - —■ r • (D8) 

xT  -C  3XJ 

It  is  worth  noting  though  that  for  x -*■  0 and  — = e <<  1 

XI 

an  asymptotic  expansion,  to  third  order,  of  the  term 
Jln(l-e)  in  equation  (D7)  identically  reduces  equation 
(D7)  to  equation  (D8) . 

In  practical  problems,  k must  be  determined  by 
patching  at  every  iteration  the  numerical  solution  around 
Xj  with  the  asymptotic  analytical  expression  (D3)  for  f(£). 
A numerical  account  of  expressions  (D7)  and  (D8)  is 
nevertheless  extremely  straightforward  if  xT  is  fixed. 

The  total  right  hand  side  of  equations  (D7)  and  (D8) 
is  evaluated  once  and  for  all  at  every  grid  point  x^ 
and  stored  in  a vector  for  successive  use. 
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APPENDIX  E 

CAUCHY  INTEGRAL  OF  A DISCONTINUOUS  FUNCTION 

Given  a function  f(x),  discontinuous  across  the  point 
x = 0 (for  convenience)  the  value  of  the  following 
integral 

a 


CI(x)  = / 

-a 


f (g)dg 
x-C 


(El) 


In  order  to  further  verify  the  result  (E2)  and  at 


the  same  time  eliminate  the  indeterminacy  of  f(x)  at 
x = 0,  a second  function  f(x)  is  also  considered,  which 
is  constant  from  -a  to  -b,  varies  linearly  from  R to  S in 
the  interval  -b  to  b and  is  again  constant  from  b to  a 
(see  figure  E2) 


Figure  E2 . Linearly  Varying  Function. 


For  this  case  1(0)  = <£  can  be  easily  determined 

-a 

analytically  and  its  limit  (for  b -*•  0)  will  produce  the 
sought  after  discontinuous  function  result,  as  follows: 

r O IS-K;  K+£>ljr.  , 3 

Jx  <f  r-TZ Tr-JdS  + 

5 


a r ~b  b (S-R)  R+S 

cko)  - $ , i isf?  + t t-2b 


-Sdi 


-a 


-a 


-b 


25  1«  + l =*? 


(E3) 
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lim  CKO)  = lim  (R-S)An  ^ = lim  (R-S)  in  . (E4) 

b+O  b+O  e-*-0  e 

If  the  function  f(x)  is  allowed  to  vary  linearly  in 
the  ranges  -a,  -b,  and  b,  a the  basic  result  does  not 
change  and  the  Cauchy  integral  CI(x)  becomes  infinite 
as  x -*•  0 like  Af  £n|x|,  (where  Af,  = S - R,  is  the  jump 
in  f(x)).  If  CI(x)  is  numerically  evaluated  by  the  second 
order  Cauchy  integral  solver  of  Appendix  A,  the  interval 
-b,  b being  equal  to  a step  size  Ax,  equation  (B6)  will 
provide  the  exact  analytical  value  of  CI(x)  at  any  grid 
point  x^  and  for  any  Ax  for  a piecewise  linear  function 
f(x).  Furthermore,  since  such  an  exact  result  is  given 
in  terms  of 
of  equation 
limit  value 
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the  same  logarithms  as  the  analytical  one 
(E3) , it  will  also  obviously  produce  the  same 
Af  in  Ax  as  Ax  0. 


L 
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Supersonic  Flow  Past  a Parabolic  Hump;  Normalized  Displacement 
Thickness  Dependence  on  h. 


85 


DDC 


-SMITH'S  LINEAR  SOLUTION 


Subsonic  Flow  Past  a Parabolic  Hump;  Normalized  Wall  Shear 


SMITH'S  LINEAR  SOLUTION 


Fig.  11.  Subsonic  Flow  Past  a Parabolic  Hump;  Normalized  Pressure 
Profile  for  h = 0.1. 
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Fig.  12.  Subsonic  Flow  Past  a Parabolic  Hump;  Longitudinal 
Step  Size  Study  for  Wall  Shear. 
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Fig.  24.  Cauchy  Integral  Model  Problem  Step  Size  Study 
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